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1.  INTRODUCTION 


The  primary  reference  surface  for  heights  is  the  geoid.  The  geoid 
is  defined  as  the  equipotential  surface  of  earth's  actual  gravity  field 
which  most  closely  approximates  mean  sea  level  in  the  spatial  least 
squares  sense.  Historically  the  geoid  has  played  an  important  conceptual 
role  in  vertical  datum  work.  A  practical  realization  of  the  geoid 
was  brought  about  by  tide  gauge  measurements  at  selected  points.  The 
effects  of  tides,  and  of  other  periodic  variations,  can  be  filtered 
out  using  appropriate  numerical  filters.  The  result  is  supposed  to 
yield  the  mean  sea  level  at  a  given  site.  Now,  to  establish  a  vertical 
network  origin  it  was  only  necessary  to  define  this  point  to  have  an 
elevation  zero.  From  this  point  levelling  can  proceed  to  determine 
potential  differences,  or  orthometric  heights,  of  all  other  points 
in  the  vertical  datum.  The  heights  then  refer  to  a  specific  vertical 
datum  by  the  mean  sea  level  at  the  origin. 

The  usual  assumption  made  in  the  past  was  that  the  mean  sea  level 
was  a  unique  equipotential  surface.  It  is  now  well  recognized  that 
this  is  not  the  case  and  that  the  vertical  datums  of  the  world  are 
not  in  a  uniform  system  due  to  the  effects  of  sea  surface  topography. 

The  sea  surface  topography  is  the  departure  of  the  mean  sea  level  from 
the  geoid.  The  causes  of  sea  surface  topography  include  ocean  currents, 
water  density  variations  as  well  as  air  pressure  and  wind  stress. 

In  addition,  close  to  the  shore  the  sea-bed  topography  and  river  discharge 
may  also  play  a  significant  role. 

It  thus  follows  that  heights  given  on  one  vertical  datum  are  not 
compatible  with  heights  on  another  datum.  Fortunately,  this  incompatibility 
is  not  large  since  the  effects  of  sea  surface  topography  ate  of  the 
order  of  2  meters.  Nevertheless,  there  are  applications  where  a  consistent 
height  system  to  a  decimeter  level  is  desired.  For  example,  in  the 
calculation  of  free  air  gravity  anomalies  to  match  today's  accuracy 
of  gravity  measurments  (0.02  mgal)  we  should  know  the  height  of  the 
point  to  6  cm  which  is  a  much  stronger  requirement  than  can  currently 
be  reached  due  to  the  variety  of  local  vertical  datums. 

The  variable  sea  surface  topography  can  cause  problems  even  within 
a  single  vertical  datum.  It  has  been  the  practice  in  some  countries 
to  connect  precise  levelling  nets  to  several  tide  gauges,  widely  dispersed 
around  the  perimeter  of  the  net.  The  sea  surface  topography  at  each 
tide  gauge  is  different,  and  resulting  adjustment  of  the  levelling 
data  to  fit  these  different  mean  sea  levels  can  lead,  and  led,  to  dis¬ 
tortions  in  the  levelling  networks. 

In  many  computations  a  precise  knowledge  of  a  set  of  spherical 
harmonic  coefficients  of  the  earth's  disturbing  potential  T  is  required. 
The  coefficients  can  be  derived  from  a  combination  of  satellite  derived 
potential  coefficients,  l°x  1°  terrestrial  gravity  anomalies,  and  l°x  1° 
anomalies  derived  from  Seasat  altimeter  data.  In  the  case  of  the  terrestrial 
data,  systematic  errors  may  be  caused  by  inconsistencies  between 


different  vertical  datums,  so  that  gravity  anomalies  after  reduction 
do  not  refer  to  a  single  equipotential  surface. 

The  purpose  of  the  present  work  is  to  evaluate  the  magnitude  of 
the  errors  in  the  potential  coefficients  caused  by  the  inconsistencies 
in  the  world's  vertical  datum.  Due  to  the  insufficient  information 
on  the  levelling  networks  in  most  of  the  areas  of  the  earth's  land 
surface,  the  technique  of  numerical  simulation  was  applied.  It  should 
be  mentioned  however,  that  the  model  height-error  function  (used  as 
the  basis  for  this  study)  was  constructed  in  such  a  way  as  to  get  the 
best  approximation  of  the  reality  from  the  information  available. 

The  contribution  of  the  computed  error  in  the  potential  coefficients 
to  the  geoid  undulations  and  gravity  anomalies  is  also  given  here. 

The  technique  of  numerical  simulation  was  also  applied  to  the 
problem  of  modelling  distortions  in  the  levelling  net  which  were  introduced 
during  the  least  squares  adjustment,  when  more  than  one  coastal  height 
was  constrianed  to  zero.  The  example  of  the  levelling  net  of  the  United 
States  was  chosen  here  because  such  procedure  actually  took  place  during 
the  1929  General  Adjustment,  and  all  heights  being  used  today  are  possibly 
affected  by  the  error  of  this  kind.  The  distortions  were  modelled 
using  the  adjustment  of  the  simulated  network  of  the  United  States 
and  were  included  in  the  basic  height-error  function. 

Summarizing,  the  results  obtained  in  these  studies  refer  to  two 
different  sources  of  errors.  The  first  is  the  error  caused  by  the  world's 
inconsistent  height  datums,  the  second  is  the  error  introduced  by  improper 
data  reduction  (forcing  the  height  of  more  than  one  coastal  point  to 
zero).  Both  sources  of  errors  are  present  in  the  actual  vertical  datum 
work,  and  both  propagate  on  to  any  derived  gravimetric  quantity  related 
to  gravity  anomalies,  such  as  the  disturbing  potential,  geoid  undulations, 
deflections  of  vertical  and  others. 


2.  SEA  SURFACE  TOPOGRAPHY  PROBLEM 

Classically  the  undisturbed  ocean  surface  has  been  taken  as  the 
geoid,  that  is  the  unique  equipotential  surface  of  the  actual  gravity 
field.  Unfortunately  such  a  surface  does  not  exist  and  the  actual 
ocean  surface  can  serve  only  as  an  approximation  of  the  geoid.  The 
departure  of  the  mean  sea  level  (MSL)  from  the  equipotential  surface 
is  called  sea  surface  topography  (SST).  The  MSL  is  constantly  changing 
with  time,  as  is  the  SST  (Merry  and  Vanicek,  1982).  The  MSL  variations 
are  mainly  due  to  salinity  differences,  temperature  differences,  atmospheri 
pressure  changes,  winds,  crustal  movements,  eustatic  changes,  and  river 
discharge.  Also,  the  global  mean  sea  surface  appears  to  be  rising 
at  a  rate  of  approximately  0.5  cm/year  (Frasetto,  1970).  Nevertheless, 
we  make  an  assumption  here  that  all  time  dependent  components  have 
been  filtered  out  and  that  SST  implied  by  such  idealized  MSL  is  independent 
of  time. 
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It  is  general  practice  to  choose  the  average  physical  water  surface 
as  the  reference  level  not  only  for  measurements  of  altitude  on  land, 
but  also  for  all  depth  determinations  at  sea.  Oceanography  studies 
of  the  MSI  require  the  determination  of  the  relative  topography  of  the 
physical  water  surface  (Lisitzin,  1974).  Different  authors  have  chosen 
a  definite  constant  depth  as  reference.  .r.  important  names  worth 
mentioning  here  are:  Dietrich  (1937),  Defant  (1941),  Lacombe  (1951), 
Reid  (1961),  Stommel  (1965),  Lisitzin  (1965),  Sturges  (1974),  and 
Levitus  (1982).  Other  oceanographers  have  chosen  a  rather  different 
approach  to  the  problem,  endeavoring  to  determine  for  every  part  of 
the  ocean  the  layer  of  'no  motion'  or,  more  precisely,  the  layer  of 
minimum  water  motion.  The  first  important  contribution  to  this  problem 
was  done  by  Dietrich  in  1937. 

Lisitzin  (1965)  computed  the  MSL  in  dynamic  cm,  relative  to  the 
4000-decibar  dynamic  depth.  According  to  the  standard  textbooks,  the 
dynamic  depth  D  of  an  isobaric  surface  is  computed  in  two  parts: 
fixed,  'standard  ocean'  part  D  (35,  0,  p)  of  uniform  35%  salinity  and 
0°  C  temperature  (for  a  given  pressure  level  p),  and  the  'anomaly  of 
dynamic  depth'  AD.  Therefore,  the  oceanic  graphics  (like  those  given 
by  Lisitzin,  1965)  represent  the  height  surplus  of  the  actual  mean 
sea  level  (or  zero-decibar  surface)  over  the  theoretical  zero-decibar 
surface  of  the  standard  ocean.  This  height  surplus  is  usually  expressed 
in  dynamic  cm,  where  (according  to  the  standard  definition)  1  dyn.cm 
is  the  unit  of  potential  corresponding  to  the  work  required  to  move 
a  unit  mass  approximately  1  cm  in  vertical  direction.  The  surface 
of  equal  dynamic  depth  is  thus  the  equi potential  surface,  although 
without  any  geodetic  significance.  Figure  1  shows  Lisitzin 's  world 
oceanic  chart  of  the  Sea  Surface  Topography.  The  sea  surface  variation 
(in  dyn.  cm)  was  computed  relative  to  the  4000-dbar  dynamic  depth  at 
open  deep  regions  of  oceans.  The  chart  reveals  that  for  all  latitudes 
the  mean  sea  level  is  lowest  in  the  Atlantic  Ocean,  while  the  highest 
values,  generally,  are  to  be  found  in  the  Pacific  Ocean.  As  an  average, 
the  Pacific  Ocean  stands  72  dyn.cm  higher  than  the  Atlantic  Ocean  and 
36  dyn.cm  higher  than  the  Indian  Ocean.  This  is  strongly  connected 
to  the  variation  of  water  density  between  the  different  oceans.  The 
highest  mean  sea  level  occurs  in  the  western  parts  of  all  oceans,  at 
the  latitudes  of  approximately  20°to  30° in  both  hemispheres.  From  this 
highest  position  the  MSL  slopes  steeply  towards  the  poles.  The  lowest 
values  occur  in  the  higher  latitudes  of  the  two  hemispheres  and  there 
is  a  secondary  minimum  in  the  zones  around  the  equator. 

The  average  sea  level  for  all  the  deep-sea  regions  of  the  oceans 
considered  in  Figure  1  is  289  dyn.cm  (Lisitzin,  1965).  The  MSL  in 
the  adjacent  and  Kaditerranean  seas  is  generally  lower  than  that  in 
the  oceans  (Lisitzin,  1965).  Therefore,  as  a  first  approximation  of 
global  average  Lisitzin  proposed  the  value  of  280  dyn.cm.  This  value 
will  be  used  in  the  sequel. 


If  we  had  reliable  values  of  SST  we  could  correct  MSL  heights 
to  get  the  geoid.  Unfortunately,  the  complexity  of  SST  problem  does 
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Figure  1.  The  distribution  of  Sea  Surface  Topography  (in  dyn.cm) 

in  the  world  oceans,  with  the  average  value  of  280  dyn.cm 
From:  (Lisitzin,  1974) 


not  allow  for  realistic  approximation  of  the  geoid  at  the  sub-decimeter 
level,  and  the  unification  of  all  vertical  datums  is  not  possible  today. 
Instead  based  on  the  existing  incomplete  information  on  SST,  we  can 
try  to  model  the  spatial  interrelations  between  different  levelling 
nets.  Extrapolating  SST  from  deep-ocean  regions  to  continental  shore-lines 
and  utilizing  available  information  on  the  levelling  techniques  and 
adjustment  procedures  « «.  different  vertical  networks  we  can  construct 
the  spatial  height-error  function.  This  function  would  reflect  the 
relative  location  of  different  vertical  datums  with  respect  to  the 
idealized  geoid.  For  the  purpose  of  this  study  the  reference  geoid 
is  that  implied  by  the  Lisitzin's  global  average  value  of  SST,  that 
is  280  dyn.cm.  The  magnitude  of  our  height-error  function  would  not 
exceed  the  amplitude  of  the  SST,  that  is  2m.  In  the  next  chapter 
we  examine  how  the  error  of  this  magnitude  propagates  on  to  the  spherical 
harmonic  coefficients  of  the  disturbing  potential  derived  from  the 
terrestrial  gravity  anomaly  data. 


3.  SST  AND  DIFFERENT  VERTICAL  NETWORK  PROCEDURES 

In  this  chapter  vertical  network  procedures  in  different  parts 
of  the  world  are  examined.  The  information  given  here  is  by  no  means 
complete.  Of  main  interest  to  us  is  the  incorrect  assumption  that 
MSL  coincides  with  the  geoid,  and  its  implication  on  the  magnitude 
of  the  relative  spatial  inconsistencies  between  different  levelling 
datums.  Such  an  estimation  was  possible  after  extrapolating  the  values 
of  SST  from  the  deep-ocean  regions  to  the  continental  shore-lines. 

This  sort  of  procedure  is  by  no  means  accurate  but  is  expected  to  provide 
a  sufficient  approximation  for  the  purpose  of  our  model  studies.  The 
Sea  Surface  Topography  on  continental  shelves  is  affected  by  the  sea-bed 
topography,  river  discharge,  wind  stress  and  other  factors  (Lisitzin, 

1965)  and,  therefore,  might  significantly  differ  from  the  values  estimated 
on  the  basis  of  extrapolation  from  the  deep  ocean  regions.  Nevertheless, 
we  use  this  technique  since  we  believe  it  assures  sufficient  accuracy 
for  the  purpose  of  this  investigation. 

The  main  source  of  information  on  SST  is  the  Lisitzin's  world-chart 
of  mean  sea  level  variations  relative  to  the  4000-dbar  isobaric  surface 
shown  on  Figure  1.  To  relate  the  extrapolated  SST  values  with  the 
vertical  datum's  zero-level  and  levelling  network  procedures,  it  is 
essential  to  know  the  spatial  distribution  of  tide-gauge  stations  which 
were  set  to  zero-height  during  the  adjustment  of  particular  levelling 
network.  Unfortunately,  such  Information  was  not  available  on  all 
continents  so  that  at  some  of  them  we  had  to  assume  that  a  unified 
continental  vertical  network  was  adjusted  with  a  hypothetical  single 
master  tide-gauge  station  constrained  to  zero-height  at  an  arbitrary 
fixed  location.  From  this  point  of  view  our  studies  are  purely  modelling 
ones  with  the  intention,  however,  to  produce  a  realistic  range  of  the 
resultant  contribution  to  the  set  of  spherical  harmonics,  and  other 
gravimetric  quantities  derived  from  gravity  data  referred  to  these 
datums.  On  the  other  hand,  wherever  possible  we  reconstruct  the  spatial 
interrelations  between  different  vertical  datums  as  realistically  as 
the  available  data  allow. 
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3.1  Classical  Vertical  Datum  Procedures  in  the  United  States  and  Canada 

The  present  vertical  reference  system  of  the  United  States  is 
based  on  the  1929  General  Adjustment  (Lippold,  1980).  This  adjustment 
contained  a  total  of  107000  km  of  levelling  in  the  U.S.  and  Canada 
(Figure  2).  The  1929  datum  was  obtained  by  holding  the  fixed  zero-elevations 
at  26  tide  gauge  sites;  21  of  these  spread  along  East,  West  and  Gulf 
Coast  of  the  U.S.,  and  5  lie  in  Canada.  The  distribution  of  primary 
tide  gauge  stations  used  as  the  reference  in  1929  adjustment  is  shown 
on  Figure  2.  The  resulting  height  system  was  designated  the  Sea  Level 
Datum  of  1929.  This  name  was  later  changed  to  the  National  Geodetic 
Vertical  Datum  of  1929  (NGVD  29). 

Canada  did  not  adopt  this  datum,  because  values  from  a  1928  adjust¬ 
ment  had  been  published  the  previous  year.  Therefore,  the  present 
Canadian  vertical  system  is  based  on  the  first  continental  adjustment 
carried  out  in  1928  and  realized  through  the  fixed  zero-height  values 
at  3  Atlantic  and  2  Pacific  tidal  stations  (Lachapelle,  1978).  Those 
five  primary  stations  plus  Rouses  Point  (located  south  of  Montreal 
on  the  U.S. -Canadian  border  and  connected  to  the  Atlantic  Ocean)  form 
the  set  of  6  zero-height  constraints  applied  during  the  1928  adjustment. 
Figure  3  shows  the  Canadian  vertical  framework  of  1928  (Lachapelle, 

1978). 

In  the  early  50's  a  new  adjustment  of  the  entire  Canadian  levelling 
network  was  made  by  the  Geodetic  Survey.  Differences  with  the  1928 
adjustment  attained  about  15  cm.,  but  no  changes  were  introduced  to 
the  1928  datum.  Summarizing,  Canadian  and  U.S.  vertical  networks  are 
on  different  vertical  reference  systems.  The  maximum  differences  reach 
10  cm  for  points  along  the  border. 

Generally,  the  latest  adjustment  of  the  North  American  Network 
in  1929  was  based  on  the  four  incorrect  assumptions  (Vanicek,  Castle, 
Balazs,  1980): 


1.  MSL  at  the  26  utilized  tide  stations  coincides  with  the  geoid 
(the  height  set  equal  to  zero); 

2.  Observed  elevation  differences  between  benchmark  pairs  are  affected 
by  only  random  error  characterized  by  a  symmetrical  distribution 

of  probability; 

3.  There  was  crustal  stability  during  the  development  of  the  network; 

4.  Orthometric  correction  computed  for  normal  gravity  instead  of 
the  actual  gravity  were  used. 

In  this  study  we  focus  only  on  the  first  assumption:  the  MSL  at  primary 
tide  stations  is  affected  by  the  irregular  SST.  Forcing  the  heights 
at  all  26  stations  to  be  zero  causes  the  distortions  in  the  resultant 
vertical  datum.  Later  we  will  evaluate  the  possible  magnitude  of  this 
distortion.  Of  course  the  distortion  will  contribute  to  an  error  in 
spherical  harmonic  coefficients  and  therefore,  will  be  incoporated 
into  our  global  model  of  height-error  function  introduced  in  Chapter  2. 

The  distribution  of  SST  values  around  the  North  American  Continent 
together  with  the  location  of  26  primary  tide  gauges  used  as  a  reference 
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Figure  2.  Primary  tide  gauge  stations  and  vertical  control  used 

in  1929  adjustment  of  the  levelling  net  in  the  United  States 
From:  (Lippold,  1980) 


Figure  3.  Canadian  vertical  framework  used  i 
the  1928  adjustment. 

From:  (Lachapelle,  1978) 


In  1929  adjustment  is  shown  in  Figure  4.  The  SST  was  extrapolated 
from  the  deep-ocean  regions  as  presented  by  Lisitzin  (1965).  It  is 
based  on  4000  dbar  isobaric  deep-sea  level.  The  units  are  dynamic 
centimeters.  The  global  average  of  280  dyn.cm  was  subtracted  from 
Lisitzin's  values  in  order  to  produce  results  with  respect  to  our  model 
geoid  implied  by  this  value. 

3.2  Vertical  Datum  Procedures  in  Australia 

On  the  Australian  continent  the  Division  of  National  Mapping  has 
conducted  a  programme  of  tidal  readings  at  30  tide  gauges  around  the 
continent  between  1966  and  1968  (Leppert,  1970).  In  1971  the  simultaneous 
adjustment  was  done  holding  MSL  for  1966-1968  fixed  at  zero  at  the 
above  mentioned  30  stations.  Figure  5  shows  the  levelling  network 
used  in  1971  adjustment  together  with  30  primary  tide  gauges  which 
were  constrained  to  zero-height.  The  net  is  split  into  5  regional 
nets  which  were  adjusted  separately. 

As  we  see,  in  the  case  of  Australian  continent  the  coincidence 
of  MSL  with  the  geoid  was  assumed.  This  incorrect  assumption  will 
lead  to  distortions  in  the  adjusted  vertical  datum.  The  distortions 
will  again  propagate  on  to  the  spherical  harmonics  coefficients.  Another 
incorrect  assumption  was  made,  when  (in  the  orthometric  correction 
formula)  the  theoretical  gravity  was  used  instead  of  the  actual  gravity. 

On  the  other  hand,  the  adjustment  programs  used  in  the  1971  procedure 
allowed  the  difference  in  height  (at  each  tide  gauge)  between  the  bench 
mark  and  sea  level  to  be  allocated  any  weight  between  zero  and  infinity. 

By  setting  weights  to  zero  Mitchell  (1972)  calculated  the  free  network 
with  a  single  constraint  at  the  fundamental  tide  station,  Jervis  Bay 
(Mitchell,  1972),  (Angus-Leppan,  Rizos,  1980).  As  will  be  pointed 
out  later,  the  only  acceptable  way  to  adjust  a  levelling  net  is  to 
use  a  single  MSL-value  constrained  to  zero-height.  Such  treatment 
leads  to  the  so  called  free  network,  which  could  be  further  block-shifted 
up  or  down  to  obtain  the  best  weighted  least  squares  fit  to  the  MSL 
at  tide  guages.  For  the  purpose  of  this  study,  we  assume  that  one 
of  the  free-network  solutions  will  be  chosen  in  the  future  to  represent 
the  vertical  datum  of  Australia.  The  extrapolated  value  of  SST  at 
the  fundamental  tide  station  in  Jervis  Bay  is  0.3  ±0.2  m  (Angus-Leppan, 
Rizos,  1980). 

Figure  6  shows  the  distribution  of  SST  values  around  the  Australian 
continent  together  with  the  location  of  30  primary  tide-station  held 
fixed  to  zero-height  during  the  1971  adustment.  The  SST  values  are 
estimated  from  Lisitzin's  map  (Lisitzin,  1965).  They  are  based  on 
the  4000  dbar  isobaric  level  and  the  global  average  of  280  dyn.cm  was 
taken  out  in  order  to  get  the  dynamic  heights  with  respect  to  our  model 
geoid  implied  by  this  value. 

3.3  Vertical  Datum  Procedure  in  Europe 


Variation' of  Sea  Surface  Topography  around 
North  America  with  respect  to  the  280  dyn.cm 
average  given  in  (Lisitzin,  1974) 

Units  are  dyn.cm. 

•-  tide  gauge  stations  used  as  a  reference 


Figure  5:  Australian  vertical  framework  of  1971 
showing  division  Into  5  regional  nets 
and  the  location  of  tide  gauges  used 
as  a  reference. 

From:  (Leppert,  1974) 


Variation  of  Sea  Surface  Topography  (Lisitzin,  1965) 
around  Australia  in  dyn.cm  (relative  to  280  dyn.cm  average) 
•  -  tide  gauge  stations 


levels  and  the  national  levelling  networks  of  Central  Europe  were  for 
the  first  time  related  to  one  common  levelling  datum  (Heer  and  Niemeier, 
1982).  In  1955,  the  International  Commission  for  the  Adjustment  of 
the  United  European  Levelling  Network  (UELN)  was  established  in  Florence. 
The  final  results  of  this  so-called  UELN-55  adjustment,  including 
Fennoscandia,  Central-  and  South-Western  Europe,  were  published  in 
1960. 

In  1981  a  combined  adjustment  of  levelling  data  from  the  continental 
countries  and  the  United  Kingdom  was  carried  out.  As  the  new  IAG  sub¬ 
commission  for  UELN  was  established  in  1973,  this  new  network  is  called 
UELN-73.  The  configuration  of  the  national  blocks  and  tidal  stations 
are  shown  on  Figure  7  (Kok,  F.,  Ehrnsperger,  W.,  Rietveld,  H.,  1980). 

To  avoid  difficulties  with  the  various  definitions  of  heights,  in  the 
adjustment  of  UELN-73  only  geopotential  numbers  or  geopotential  units 
(gpu)  were  introduced  as  parameters  representing  the  'heights’  of  all 
surface  points,  and  geopotential  differences  were  introduced  as  obser¬ 
vations.  In  phase  I  of  this  new  computation  the  free  adjustment  of 
the  levelling  networks  was  performed,  without  incorporating  any  infor¬ 
mation  on  MSL.  At  first,  the  data  of  the  14  West  European  countries 
participating  in  the  net  were  adjusted  separately.  After  these  prelim¬ 
inary  computations  a  combined  adjustment  was  carried  out,  from  which 
the  geopotential  numbers  of  the  included  reference  points  and  their 
standard  deviations  were  taken.  A  single  datum  point  in  this  adjustment 
is  the  Normal  Amsterdam  Piel  (NAP)  represented  in  UELN-73,  by  point 
No.  4019  (Kok,  Ehrnsperger,  Rietveld,  1980).  Gravity  data  was  based 
on  the  International  Gravity  Standardization  Net  1971  (IGSN  71).  The 
method  of  observation  equations  was  used.  The  79  tidal  stations  are 
connected  to  the  net,  but  the  information  on  MSL  was  not  used  as  con¬ 
straints,  allowing  for  the  free  adjustment  dependent  only  on  the  geo¬ 
potential  number  of  the  datum  point  at  NAP. 

In  phase  II  of  UELN-73  it  is  planned  to  perform  statistical  testing 
of  adjusted  heights  of  tidal  stations  against  the  oceanographic  definition 
of  MSL.  In  the  future,  based  on  such  tests  the  whole  block  of  adjusted 
values  obtained  in  the  free  adjustment  of  phase  I  could  be  shifted 
in  vertical  direction  in  order  to  obtain  an  optimal  least  squares 
fit  to  the  oceanographic  MSL-values  at  tide  stations.  After  removing 
the  average  value  of  SST  from  coastal  heights  such  optimal  fit  would 
ultimately  produce  the  geoid,  at  the  resolution  of  the  magnitude  of 
uncertainty  to  which  the  values  of  SST  are  known. 

Figure  8  shows  the  distribution  of  SST  values  around  Europe  together 
with  the  location  of  tide  gauges  connected  (passively)  to  the  net  (Kok,  J. 
Ehrnsperger,  W.  Rietveld,  H.,  1980).  The  SST  heights  are  estimated 
from  Lisitzin's  map  and  the  global  average  of  280  dyn.cm  was  taken 
out  in  order  to  produce  the  dynamic  heights  with  respect  to  our  model 
geoid  implied  by  this  value. 

Since  no  constraints  have  been  imposed  on  heights  at  coastal  point, 
the  UELN-73  datum  is  f^ee  of  internal  distortions  caused  by  the  departure 


of  MSL  from  the  geoid.  However,  the  systematic  error  exists  in  all 
values  and  is  equal  to  the  SST  value  at  the  datum  point  (NAP).  This 
value  could  be  estimated  from  Figure  8  as  about  -50  cm  (Lisitzin,  1974, 
p.  1963). 

3.4  The  SST  Estimates  Alona  Other  Countries 


The  Lisitzin's  world's  oceanic  chart  (Figure  1)  can  serve  as  a 
basis  for  extrapolation  of  the  SST  values  along  the  coasts  of  other 
continental  blocks.  The  values  of  SST  are  based  on  the  4000  dbar  iso 
baric  level  and  the  global  average  of  280  dyn.cm  was  taken  out. 

Figure  9  shows  the  distribution  of  SST  around  Asia,  Figure  10 
around  the  coasts  of  South  America,  and  Figure  11  around  Africa.  The 
units  for  Sea  Surface  Topography  variation  are  dynamic  centimeters. 


4.  THE  DISTORTIONS  IN  A  VERTICAL  DATUM  DUE  TO  THE  ZERO-CONSTRAINTS 


AT  SEVERAL  TIDE  GAUGES 

As  we  know  today,  the  classical  assumption  that  the  Mean  Sea  Level 
coincides  with  the  geoid  is  not  valid.  However,  even  today  it  remains 
a  common  practice  to  assign  the  zero-heights  to  coastal  points  connected 
to  the  levelling  net.  These  so-called  zero-constraints  introduce  inherent 
distortions  in  the  adjusted  vertical  datum.  The  magnitude  of  this 
distortion  depends  mainly  on  the  configuration  of  tide  gauges  which 
were  constrained  to  zero-height  with  respect  to  the  particular  dis¬ 
tribution  of  SST  along  the  perimeter  of  the  net.  The  distortions  intro¬ 
duced  at  tidal  stations  can,  therefore,  reach  the  amplitude  of  variation 
of  SST  along  the  ocean  boundaries  of  the  given  network.  During  the 
adjustment  process  those  distortions  propagate  inside  the  net  and  all 
the  interior  nodes  will  be  affected  by  this  displacement.  In  this 
chapter  we  want  to  evaluate  the  magnitude  of  this  displacement  on  the 
example  of  the  levelling  net  of  the  United  States. 

All  heights  on  the  North  American  continent  are  based  on  the 
National  Geodetic  Vertical  Datum  of  1929.  Therefore,  the  North  Ameri¬ 
can  vertical  datum  is  basically  warped  due  to  the  26  zero-constraints 
at  tide  gauges  incorporated  to  the  net  during  adjustment.  The  magnitude 
and  spatial  behaviour  of  this  distortion  is  of  great  importance  to 
our  present  problem.  We  will  model  this  distortion  and  next  evaluate 
its  effect  on  potential  coefficient  together  with  the  effects  resulting 
from  the  SST  problem. 

4.1  The  General  Approach  to  the  Problem 

In  order  to  evaluate  the  distortions  in  the  levelling  net  the 
method  of  numerical  simulation  was  used. 

First,  we  design  a  simplified  regular  grid  to  model  the  1929- 
levelling  net  of  the  U.S.  Our  numerical  approximations  to  the  actual 
network  will  be  constructed  as  a  regular  grid  inside  the  rectangle 
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1920  km  NS  and  4000  km  EW,  which  more  or  less  approximates  the  actual 
boundaries  of  the  U.S.  Figure  12  shows  the  example  of  such  network 
for  the  case  of  36(*6x6)  grid  points.  The  net  basically  consists  of 
25  rectangular  loops,  each  384  km  in  North-South  direction  and  800  km 
in  East-West.  Each  grid-point  represents  the  unknown  height,  therefore, 
there  are  36  unknowns  in  this  example.  In  the  actual  vertical  datum 
work  the  observations  are  taken  as  height  differences  between  points 
in  the  net.  In  the  example  we  assume  there  is  one  observation  (height 
difference)  between  each  pair  of  adjacent  points  giving  60  observation 
equations.  It  means  no  observation  stretches  over  more  than  one  segment 
in  the  net. 

As  in  the  case  of  actual  adjustment  of  the  1929  levelling  net, 
each  segment  between  'bench  marks'  in  the  example  provides  one  linear 
observation  equation  of  the  type 


or 


(1) 


Ah  =  Ah 


(2) 


In  order  to  simulate  the  set  of  observations  Ah-jj  we  first  assign 
heights  (from  a  topographic  map)  to  all  grid  points  h°,  h°,  ....  hfj. 
The  grid  points  are  ordered  as  shown  on  Figure  12,  which  assures  a 
regular  structure  of  the  linear  system  of  observation  equations.  At 
this  stage  we  assign  the  sea  surface  topography  heights  to  'coastal' 
points  marked  on  Figure  12.  The  values  of  SST  with  respect  to 
280  dyn.cm  global  average  can  be  estimated  from  Lisitzin's  chart 
(Figure  4).  Now,  starting  from  those  reference  heights,  we  can  compute 
the  set  of  observations  Ahjj  using  (1).  At  this  point  we  have  the 
right-hand  side  and  the  design  matrix  A  of  formula  (2).  The  standard 
deviations  of  our  synthetic  observations  were  assigned  according  to 
the  usual  formula  (Milbert,  1980): 


where  Li j  -  section  length  (in  km)  between  bench  marks  i  and  j 
and  oQ  =  0.3  cm.  In  our  simple  example  it  gives  the  values: 

8.49  cm  in  East-West  direction 

°ij  '  ' 

5.88  cm  in  North-South  direction 

V 
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Now,  we  solve  the  system  of  observation  equations  (2)  using  the 
method  of  linear  least  squares  with  linear  equality  constraints  by 
downweighting.  A  brief  description  of  the  method  together  with  some 
accuracy  considerations  will  be  given  in  the  next  section.  We  can 
obtain  the  solution  with  different  sets  of  constraints.  Since  in  the 
1929  adjustment  of  the  levelling  net  of  the  U.S.  each  geodetically 
connected  tide  gauge  provided  a  zero-height  constraint,  in  our  model- 
studies  we  also  attached  a  set  of  linear  constraints  of  the  type: 


hj  “0,  i  =  1 ,  . . . ,  k 


to  the  grid  points  representing  ‘tide  gauges’  along  the  West  Coast, 
the  East  Coast,  and  the  Gulf  Coast.  In  our  6  by  6  example  we  have 
13  constraints  equations,  forcing  the  least  squares  solution  to  assume 
zero-height  values  at  13^grid  points  marked  on  Figure  12.  Let  this 
solution  be  denoted  by  h  . 

Now  we  can  subtract  our  reference  "true"  heights  h°  from  the 
constrained  least  squares  solution  h: 

d  =  h  -  h°  (3) 


The  function  d  represents  the  vertical  datum  distortion  caused  by 
the  zero-constraints  imposed  on  heights  at  tidal  stations.  This  function 
should  be  understood  in  the  spatial  sense  -  the  distortion  varies  from 
point  to  point  forming  a  more  or  less  irregular  but  continuous  surface. 

The  shape  and  the  amplitude  of  this  function  is  of  particular  importance. 

We  can  assume  that  heights  in  the  United  States,  and  most  likely  in 
the  entire  N.  American  Continent  that  are  used  today,  are  contaminated 
by  distortion  of  this  type.  The  North  American  vertical  datum  is  warped 
and  this  deformation  can  be  quantitatively  modeled  by  our  function  d. 

The  contribution  of  this  distortion  to  the  error  in  spherical  harmonic 
coefficients  can  be  now  evaluated  on  performing  the  Fourier  analysis 
of  the  distortion  function  d  as  will  be  described  later. 

4.2  The  Software  Considerations 

A  method  of  linear  least  squares  with  linear  equality  constraints 
by  downweighting  was  used  to  solve  for  the  distortion  function  d 
in  our  model  studies.  This  method  which  is  based  on  orthogonal  decomposition 
is  numerically  more  stable  than  the  usual  methods  based  on  normal  equations 
approach.  For  this  reason  it  was  used  to  handle  our  problem.  The  House¬ 
holder  orthogonal  decomposition  plays  a  major  role  in  this  technique, 
therefore  we  start  from  a  brief  summary  of  some  algebraic  principles 
behind  this  approach. 

4.2-1  Householder  Orthogonal  Decomposition  for  Least  Squares  Problems 
and  Disadvantages  of  Normal  Equations  Method 


An  orthgonal  decomposition  of  mxn  matrix  A  of  rank  k  means 
to  find  an  mxm  orthogonal  matrix  H,  an  mxn  matrix  of  the  form 


0 

0 


and  an  nxn  orthogonal  matrix  K,  so  that  A=HRKT  and  Rn  is  a 
kxk  matrix  of  rank  k. 


An  important  property  of  the  orthogonal  matrix:  if  Q  is  ortho 
gonal,  then  we  have: 


-  preservation  of  Euclidean  length  (norm)  under  multiplication 

IIQylMM! 

-  stability  in  computation  with  regard  to  propagating  data  errors  or 
uncertainty. 


Householder  transformation  (Householder,  1958). 

Given  a  m-vector  v^O,  there  exists  an  orthogonal  matrix  Q  such  that 


Qv  =  -a  II  v ||  e  i  with 


Y 

Vf 

0 

S' 

V2 

+1  if  Vj  >  0 

and  a  =  - 

where  v  = 

• 

-1  if  v,  <  0 

* 

• 

^  l 

# 

0, 

vm 

It  turns  out  that  Q  =  I - y—  ,  where  u  =  v  +  o  ||v||  e1(  and 

Im  -  identity  matrix.  u  u 

Geometric  interpretation  of  Householder  transformation: 

it  represents  a  reflection  in  the  (m-1)  -  dimensional  subspace,  S  , 
orthogonal  to  the  vector  u  .  By  this  it  is  meant  that  Qu  =  u  and 
Qs  *  s  for  all  s  e  S. 

At  this  stage,  we  can  formulate  the  least  squares  problem  (Lawson 
and  Hanson,  1974). 


For  an  equation 


with  a  real  mxn  matrix  A  of  rank  k  min  (m,m)  and  a  real  m-vector 
b  ,  find  a  real  n-vector  x  minimizing  the  euclidean  length  of  Ax-b. 
Additionally,  we  explicitly  assume  that  the  numerical  data  that  constitute 
A  and  b  have  only  a  limited  number  of  accurate  digits.  This  assumption 
reflects  the  usual  situation  in  numerical  handling  of  data.  In  our 
modelling  studies  the  set  of  observation  equations  (2)  will  be  considered 
a  least  squares  problem. 

Next  we  briefly  describe  the  numerical  properties  of  the  algorithm 
to  solve  problem  least  squares  by  orthogonal  decomposition  using  House¬ 
holder  transformations  (Golub  and  Businger  (1965),  Hanson  and  Lawson 
(1969)).  Depending  on  the  rank  of  the  design  matrix  A  of  a  least 
squares  problem  Ax  =  b  the  algorithm  takes  the  following  actions: 

-  if  matrix  A  is  well  conditioned  (full  rank)  this  algorithm  finds 
a  standard,  unique  least  squares  solution  x  that  minimizes 

Ax  -  b  . 

-  if  matrix  A  is  rank-deficient  the  algorithm  determines  the  effective 
rank  k  of  matrix  A  and  finds  the  so  called  minimum  length  solution 
x  that  minimizes  both:  ||Ax  -  b({  and  ||x||.  This  solution  is  also 
unique.  The  effective  rank  k  is  the  rank  of  the  matrix  A  that  re¬ 
places  A  as  a  result  of  a  specific  computational  algorithm. 

-  if  matrix  A  has  theoretically  full  rank,  but  inevitable  small  changes 
in  the  data  of  the  order  of  magnitude  of  data  uncertainty  could  convert 
the  matrix  to  one  of  deficient  rank  thecal gorithm  replaces  A  by 

a  'nearby'  rank-deficient  matrix,  say,  A  ,  (Lawson  and  Hanson,  1974, 
p.  79,  eq.  14.8),  and  then  competes  a  minimum  length  solution  to  the 
slightly  different  problem  Ax  ?  b. 

The  last  case  is  the  ill-conditioned  case  and  is  traditionally 
disregarded  by  standard  algorithms  that  solve  least  squares  problems. 

Indeed,  it  is  a  difficult  problem  to  handle  numerically,  and  (if  not 
properly  treated)  can  lead  to  very  unstable,  inaccurate  solutions. 

The  rank  of  the  matrix  A  that  replaces  A  as  a  result  of  a  specific 
computational  procedure  is  called  the  pseudorank  of  A.  Pseudorank  (or  ef¬ 
fective  rank)  is  not  a  unique  property  of  Matrix  A,  but  also  depends  on  the 
details  of  computational  algorithm,  the  values  of  assumed  tolerance  parameters 
and  the  effects  of  machine  round-off  errors. 

Generally,  the  algorithm  (Lawson  and  Hanson,  1974,  p.  78)  determines 
the  orthogonal  matrices  Q  ,  K  and  the  permutation  matrix  P  in  the 
following  steps: 

DECOMPOSITION 

a)  find  orthogonal  matrices  Q,  P  such  that 
A  *  QtR  PT,  or 


k  n-k 
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Wyi  =  d 

and  set  y2  arbitrary 

b)  the  solution  vector  x  .is  found  as  follows: 

K  0  1  Vil  )  'k 

X  *  p 

0  ln-k  *2  )  n’k 


c)  the  residual  norm  can  then  be  computed  as  ||b-Ax||  =||c2||. 


d)  the  variance-covariance  matrix  of  the  solution  vector  x  for  the 
full  rank  case  (k=n)  is  proportional  to 


C:  =  PR-1  (R-1)1  PT 


(Lawson  and  Hanson,  1974,  p.  68,  eq.  12.10) 

Both  orthogonal  matrices  Q  and  K  are  products  of  a  number  of  elementary 
Householder  transformations. 

Remarks  on  the  algorithms: 

-  The  described  algorithm  solves  directly  the  original  least  squares 
problem,  it  means  the  system  of  observation  equation  with  rectangular 
design  matrix  A:  Ax^o.  (the  formulation  of  the  so-called  normal 
equations  is  avoided). 

-  Because  the  orthogonal  decomposition  of  A  is  found,  there  is  no 
(explicit)  matrix  -  inversion. 

One  may  think  that  using  a  general  inverses  approach  (presented 
by  Householder  method)  one  has  to  resign  from  statistical  indicators 
as  covariance  matrices  (which  are  formally  related  to  the  normal  matrices 
(symmetric,  positive  definite)).  This  is  not  true,  as  can  be  found 
in  Lawson  and  Hanson  (1974,  pp.  67-73).  It  turns  out  that  using  the 
computed  components  of  the  Householder's  orthogonal  decomposition, 
one  can  construct  the  required  covariance  matrices  in  a  cheap,  effec¬ 
tive  and  stable  way. 

4.2.2  Method  of  Linear  Least  Squares  with  Linear  Equality  Constraints 
by  Downweighting. 

During  the  1929  adjustment  of  the  levelling  net  of  the  U.S.  the 
heights  at  each  of  the  26  geodetically  connected  tide  gauges  were  con¬ 
strained  to  zero  through  a  linear  equation  of  the  type 


h-j  =  0,  i  =  l,  2,  ....  k  (hi  stands  for  the  height  at  tide  station  i) 


or  generally 


Ch  =  f 


(5) 


In  the  adjustment  of  our  model -net  we  also  wish  to  impose  zero- 
constraints  on  grid  points  bordering  upon  the  sea.  As  can  be  found 
in  Lawson  and  Hanson  (1974),  this  can  be  done  simply  by  entering  the 
associated  constraints  equations  (5)  before  our  original  observation 


equation  system  (2)  or  (4).  Next,  the  set  of  the  original  observation 
equations  (4)  should  be  downweighted  by  premultiplying  by  some  small 
quantity  e  to  be  chosen  by  the  user.  This  quantity  establishes  the 
compromise  between  the  usual  least  squares  solution  of  the  unconstrained 
system  (4)  and  the  closeness  of  fit  to  the  prescribed  constraints  (5). 
However,  the  internal  relative  importance  of  each  equation  in  the  un¬ 
constrained  system  (4)  expressed  by  the  weights  associated  with  the 
original  observations  (Section  4.1)  is  preserved  during  this  down¬ 
weighting  operation.  Summarizing,  we  tombine  observation  equations 
Ah  =  Ah  and  constraints  equations  Ch  =  f  into  one  linear  system: 


c" 

f 

h  = 

^  <*> 

i _ 

e  Ah 

(6) 


This  system  can  be  solved  for  the  height  vector  h  using  the  method 
of  Householder  orthogonal  decomposition  described  in  the  previous 
section. 

It  should  be  mentioned,  that  in  our  case  of  zero-constraints  vector 
f  in  (6)  is  set  to  0. 

Altogether,  this  numerical  procedure  is  called  in  the  literature, 
the  least  squares  method  with  linear  equality  constraints  by  down¬ 
weighting  (Lawson  and  Hanson,  1974,  p.  148-158).  And  as  mentioned 
already,  the  small  parameter  e  expresses  the  compromise  between  the 
closeness  of  fit  to  the  values  prescribed  by  the  constraints  and  fi¬ 
delity  to  the  unconstrained  least  squares  solution  of  the  original 
system. 


Ah  =  Ah 


Powell  and  Reid  (1968)  have  proven  the  stability  of  this  method,  provided 
the  Householder  transformations  are  used,  and  constraint  equations 
are  placed  on  the  top  of  the  original  system  (2). 

4.2.3  Comparison  with  Normal  Equations  Method 

It  can  be  shown  (Lawson  and  Hanson  (1974),  p.  122)  that  for  a  mxn 
matrix  A  the  Householder  method  applied  directly  to  least  square  prob¬ 
lem  Ax  =  b,  requires  mn2  +  m3/3  operations  (by  an  operation  we  mean  a 
pair:  addition  -  multiplication).  The  traditional  method  by  normal 
equations  requires  mn2/2  operations  to  form  the  normal  equations 
A^Ax  =  A'b,  and  n3/6  operations  to  solve  normal  equations  using 
Cholesky's  method.  Therefore,  Householder  method  is  twice  as  expensive 
as  the  traditional  normal  equation  method.  On  the  other  hand  using 
computer  arithmetic  of  the  given  precision  the  Householder  method 
produces  twice  as  many  significant  digits  in  the  solution  as  the 
normal  equation  method  does.  In  other  words,  to  assure  the  same  quality 
of  results  we  have  to  double  the  computer  precision  while  using 
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normal  equations  method.  For  large  systems  of  equations  it  leads  quickly 
to  computer-storage  problems. 


Conclusion:  on  a  given  computer  the  Householder  method  can  handle 
larger  problems  than  the  normal  equation  method  can. 

Lawson  and  Hanson  (1974)  give  the  following  simple  example  of 
the  loss  of  significant  digits  using  normal  equations  method: 


1 

1 


1  1-e  J 


The  normal  equations  matrix  A  is  to  be  formed  on  the  computer 
with  the  relative  precision  n  .  Suppose  the  value  of  e  is  such, 
that  it  is  significant  to  the  problem,  e  >  lOOn  say,  but  e2  <  n  , 
so  that  1-e  f  1  but  3+e2  is  computed  as  3.  Thus,  instead  of  computing 


aTa  = 


3 

3-e 


3-e 

3-2e+e2 


we  shall  compute  the  matrix 

3  3-e 

3-e  3-2e 


Further  computation  of  Cholesky  decomposition  of  this  matrix  into 
the  product  R^R  ,  gives 


where 

rn  =  /J,  rl2  =  (3-e )//3,  r22  B  3-2e-  (9-6e)/3  =  0,  while  the  correct 

result  at  this  point  would  be 


r\ 2  =  2  e2/3  /  0. 
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Therefore,  in  the  precision  n  no  significant  digits  are  obtained 
in  the  element  r22  •  Consequently  the  matrix  R  is  computed  as 
singular.  On  the  other  hand,  when  the  Householder  transformation  is 
applied  to  triangul arize  A  directly,  without  forming  A^A  ,  we  get 
the  matrix 


e-3  ~l 

/I 

-€  /E 
3 


which  is  nonsingular. 

The  important  point  in  the  Householder  routine  is  that  components 
of  R  were  computed  as  the  difference  of  quantities  of  order  of  magnitude 
unity.  Thus,  using  n-precision  arithmetic,  these  components  are  not 
lost  in  the  round-off  errors. 

4.3  The  Results  of  Simulation  of  the  North  American  Vertical  Datum 
Distortion  Function  for  Different  Densities  of  Model  Networks 

In  Figure  12  we  see  the  example  of  a  synthetic  network  for  the 
array  of  6x6  grid  points.  Points  are  regularly  distributed  and  ordered 
row-wise  so  that  the  densification  of  the  net  could  be  easily  programmed. 
The  least  squares  method  with  linear  constraints  by  downweighting  was 
used  to  obtain  the  solution  to  the  combined  system  of  observation  equa¬ 
tions  together  with  constraints  equations  as  was  described  in  previous 
section.  We  denote  this  solution  by  h  ,  the  actual  dimensions  of 
linear  system  to  be  solved  were  for  this  case:  36  unknown  heights 
and  73  equations  (60  observations  +  13  constraints).  The  difference: 
constrained  solution  h  minus  the  original  reference  heights  h° 
described  in  Section  4.1,  gives  the  distortion  function  of  the  vertical 
datum. 


d  =  h  -  h° 


Figure  13  shows  the  contour  map  of  the  distortion  function  d 
for  the  6x6  network,  together  with  the  3-dimensional  view  of  this  func¬ 
tion.  The  units  are  cm  and  contour  interval  is  equal  to  2  cm.  The 
downweighting  parameter  e  has  the  value  0.1  which  is  sufficient  to 
constrain  the  zero-height  at  tide  gauges  with  tolerance  of  2  mm.  From 
the  contour  map  we  see  that  the  distortions  assume  the  extremal  values 
at  tide  stations,  where  the  zero-constraints  were  forced.  The  minimum 
occurs  at  the  south  part  of  the  West  Coast,  the  maximum  at  East  Coast, 
near  Florida.  From  South-West  Coast  the  distortion  function  increases 
gradually  toward  East.  Approximately  in  the  middle  of  the  continent 


Figure  12.  Example  of  the  6x6  rectangular  grid  approximating  the  1929  Levelling  Net 
in  the  U.S.A.,  showing  the  order  of  unknowns  and  the  constrained  heiqhts 
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the  distortions  diminish  to  zero  and  the  line  of  no  distortion  passes 
across  the  continent  almost  exactly  in  the  North-South  direction. 

From  here  distortions  increase  again  towards  East  and  South.  Of  course 
this  particular  shape  of  the  distortion  function  d  is  a  result  of 
a  specific  distribution  of  sea  surface  topography  around  North-American 
c  itinent  shown  on  Figure  4.  It  can  not  be  generalized  to  any  other 
continent  or  location. 

In  the  next  stage  we  densify  the  model  levelling-network.  For 
the  case  of  a  10x10  regular  grid  the  dimensions  of  the  linear  system 
are:  100  unknown  heights,  and  204  equations  (180  observations  +  24 
constraints).  Figure  14  shows  the  contour  map  and  the  3  dimensional 
view  of  the  distortion  function  d  for  this  case.  The  contour  inter¬ 
val  is  equal  to  2  cm  and  the  downweighting  parameter  e=0.1,  which 
is  enough  to  constrain  the  heights  to  zero  at  24  'tide  gauges'  with 
the  tolerance  of  1  mm.  The  distribution  of  the  distortions  over  the 
continent  is  very  similar  to  the  previous  case. 

Finally,  in  Figure  15  the  contour  map  and  the  3-dimensional  view 
of  the  distortion  function  for  the  15x  15  grid  is  shown.  Units  are 
cm  and  contour  interval  is  equal  to  2  cm.  The  downweighting  parameter 
e  is  equal  to  0.2  which  is  sufficient  to  constrain  the  heights  to 
zero  at  35  'tide  gauges'  with  the  tolerance  of  1  mm.  For  this  case 
the  system  of  445  equations  with  225  unknowns  was  solved.  The  general 
behaviour  of  distortion  function  is  again  similar  to  the  previous  ex¬ 
amples. 

Summarizing  the  above  results  we  can  formulate  the  following  general 
remarks: 

-  the  maximum  (in  magnitude)  distortions  in  the  vertical  datum  occur 
at  the  tide  gauges,  where  incorrect  zero-constraints  were  forced: 

-  the  extreme  distortions  values  are:  -20  cm  at  the  South  part  of 
West  Coast  and  +21  cm  at  the  South  Part  of  the  East  Coast; 

-  in  the  center  of  the  continent  the  distortions  never  exceed  those 
at  the  tide  gauges; 

-  the  shape  of  the  distortion  function  is  not  sensitive  to  the  density 
of  grid,  it  is  a  result  of  a  particular  distribution  of  the  SST  along 
the  coasts  of  the  continent; 

-  the  line  of  no  distortion  passes  through  the  middle  of  the  continent 
in  the  meridional  direction; 

-  the  vertical  datum  was  warped  in  a  smooth  way,  which  means  there 

is  no  unexpected  extremes  (local  minimums  or  maximums)  in  the  center 
of  continent,  and  there  is  no  long  or  short-wavelength  periodic  effects. 

The  last  remark  may  come  from  the  fact  that  our  synthetic  network 
was  very  sparse,  which  is  equivalent  to  the  smoothing  effect.  On  the 
other  hand,  the  general  shape  of  continental  topography  (which  is  imaged 
in  the  input  data)  is  very  unsymmetric,  with  high  mountains  on  West 
Side  and  plains  on  East.  Nevertheless,  the  topography  does  not  seem 
to  affect  the  shape  of  the  distortion  function. 


,V>  f  t  {  /  t  * 


Of  particular  importance  is  the  existence  of  the  zone  of  no  dis¬ 
tortion.  It  indicates  that  along  North-South  span,  right  in  the  middle 
of  the  continent  the  heights  are  not  affected  by  error  due  to  the 
twist  in  the  vertical  datum.  In  other  words,  the  equipotential  surface 
implied  by  the  zero-height  for  points  inside  this  zone  could  represent 
the  mean  regional  geoid  for  the  North  American  continent.  At  all  other 
points  heights  refer  to  the  different  equipotential  surfaces  which 
deviate  from  the  mean  regional  geoid  exactly  by  the  particular  value 
of  the  distortion  function  d  .  It  means  that  on  adding  the  distortion 
function  d  ,  all  heights  on  the  continent  would  refer  to  a  single 
equipotential  surface. 


4.4  The  Proper  Procedure  for  the  Adjustment  of  a  Levelling  Network 
and  the  Definition  of  the  Geoid 


In  Section  4.2  it  was  shown  that  forcing  the  zero-height  at  more 
than  one  tide  gauge  introduces  errors  in  the  adjusted  heights  of  the 
order  of  the  variation  of  SST  at  tide  stations  chosen.  In  order  to 
avoid  this  effect  we  should  allow  for  a  minimum-constraints  solution, 
where  a  single  point  is  assigned  an  a  priori  value.  After  the  network 
is  adjusted  the  computed  heights  of  tidal  stations  should  be  statis¬ 
tically  tested  against  the  MSL  heights  determined  by  the  oceanographic 
method.  Next,  the  vertical  translation  of  all  adjusted  heights  could 
be  designed  in  order  to  minimize  the  differences  between  the  MSL  values 
and  the  adjusted  heights  at  tidal  stations  in  some  least  squares  sense. 
Heights  obtained  that  way  would  refer  to  a  single  equipotential  surface. 
This  surface  could  be  regarded  as  the  definition  of  the  geoid  in  the 
regional  sense.  It  would  have  the  property  of  minimizing  the  Sea  Surface 
Topography  variations  over  the  given  set  of  tidal  stations.  Therefore, 
such  a  regional  geoid  would  depend  on  the  particular  configuration 
of  tide  gauges  used.  The  conceptual  difference  between  the  regional 
geoid  introduced  here  and  the  geoid  itself  is  the  following:  The  geoid 
(at  least  in  the  oceanographic  sense)  is  the  equipotential  surface 
of  the  earth's  gravity  field  that  most  closely  (in  some  more  or  less 
formal  i zed  sense)  approximates  the  MSL  in  the  global  context.  The  region¬ 
al  geoid  (for  a  given  continent)  is  that  equipotential  surface  of 
the  earth's  gravity  field  that  most  closely  approximates  in  the  least 
squares  sense  the  MSL  variation  over  the  discrete  set  of  coastal  points 
(tide  gauge  locations).  The  vertical  datum  implied  by  the  regional 
geoid  would  define  the  zero-level  that  on  average  agrees  with  the  MSL 
along  the  oceanic  coasts  of  the  region.  Essentially  one-dimensional 
information  on  the  MSL  variation  along  the  coastal  line  (sampled  at 
the  tide  gauge  locations)  is  needed  to  define  the  regional  geoid, 
whereas  the  areal  distribution  in  the  two-dimensional  global  sense 
is  required  to  define  geoid. 


Another  way  to  avoid  distortions  in  the  adjusted  vertical  datum 
is  to  perform  an  unconstrained  adjustment.  Of  course  the  linear  system 
that  arises  from  this  problem  is  singular  and  can  not  be  solved  by 
standard  methods  of  normal  equations.  On  the  other  hand  the  House¬ 
holder  orthogonal  decomposition  (presented  in  Section  4.2.1)  is  well 
suited  for  rank-deficient  systems.  In  such  case  the  algorithm  produces 


the  so-called  minimum  length  solution  which  is  unique  in  this  sense. 
Next,  this  solution  could  be  block-shifted  in  vertical  direction  to 
meet  certain  conditions  at  tide  gauges,  which  were  described  in  the 
last  paragraph. 

In  order  to  find  out  the  differences  between  those  two  methods 
a  synthetic  network  of  7  by  7  grid  points  was  adjusted  exactly  the 
same  way  as  described  in  Section  4.1.  The  only  difference  was,  that 
a  single  zero-constraint  was  imposed  on  point  at  the  North-West  corner 
of  the  net,  just  to  keep  the  full  rank  of  linear  system.  At  the  same 
time  an  unconstrained  minimum  length  solution  was  obtained  for  the 
same  case  of  the  7x7  network.  Then,  this  solution  was  block-shifted 
so  that  the  height  at  North-West  corner  of  the  net  became  zero 
(exactly  like  in  the  single-constraint  solution).  Figure  16  shows 
the  difference  between  the  n^-constraints  and  the  single-constraint 
solution.  The  units  are  10"  ^  m.  The  maximum  difference  was  4-10"12  m, 
which  indicates  that  solutions  differ  only  on  the  order  of  round-off 
errors. 

The  above  example  indicates  that  both  methods  produce  equivalent 
results.  Either  one  can  be  applied  in  real-life  situations. 


5.  THE  DESIGN  OF  THE  MODEL  HEIGHT-ERROR  FUNCTION  DUE  TO  THE 
INCONSISTENCIES  IN  THE  GLOBAL  VERTICAL  DAtUM 


Now  we  are  in  the  position  of  constructiong  the  global  model  of 
the  height-inconsistency  function  Ah(  4>,  x)  caused  by  the  inconsis¬ 
tency  between  the  different  vertical  datums  in  the  world.  In  the 
spatial  sense  ah(<t> ,  x)  is  understood  to  be  a  function  of  geodetic 
latitude  <t>  and  longitude  x  .  We  are  interested  only  in  the  effects 
related  to  the  uncertainties  in  the  Sea  Surface  Topography  values  or 
due  to  the  improper  numerical  procedure  during  the  adjustment  of 
levelling  networks. 


The  reference  surface  for  the  construction  of  Ah(  x)  is  the 
equipotential  surface  implied  by  the  Lisitzin's  global  average  value 
of  280  dyn.cm  of  the  SST  variation  (Lisitzin,  1965).  This  value  was 
derived  based  on  the  4000  dbar  isobaric  depth. 


The  function  Ah  is  modeled  as  the  patchwise  step  function  as¬ 
suming  constant  values  over  specific  vertical  datums.  Such  constant 
values  refer  to  the  sea  surface  topography  at  a  particular  tide  gauge 
which  was  constrained  to  zero  during  the  general  adjustment  of  level¬ 
ling  net.  In  case  where  more  than  one  tide  gauge  were  constrained 
to  zero  (N.  America)  the  average  sea  surface  topography  at  all  such 
tide  gauges  was  chosen  to  represent  the  correction  Ah  over  this  par¬ 
ticular  vertical  datum.  This  value  was  estimated  based  on  the  partic¬ 
ular  distribution  of  SST  around  the  perimeter  of  the  continental  level¬ 
ling  network.  The  Lisitzin's  deep-ocean  Sea  Surface  Topography  as 
shown  on  Figure  1  was  extrapolated  towards  the  locations  of  tide  gauges 
which  took  part  in  the  adjustment  of  the  levelling  net. 


The  difference  between  the  unconstrained 
(rank-deficient)  solution  and  the  minimum 
constraints  (full-rank)  solution  for  the 
case  of  7x7  Digital  Surface  Model. 

Units  are  pm=10"I2m and  contour  interval  is  0.3 


At  the  ocean  areas  the  Lisitzin's  Sea  Surface  Topography  values 
are  assumed  to  represent  the  discrepancy  Ah  between  the  MSL  and  our 
reference  geoid  implied  by  the  280  dyn.cm  global  average. 

More  specifically,  over  the  continental  plates  we  imposed  the 
following  constant  values  for  Ah(4>,  a): 

a)  North  America 

In  the  1929  levelling  net  adjsustment  26  Primary  Tide  Stations  were 
used  as  the  reference  zero-level.  The  average  value  of  Sea  Surface 
Topography  over  those  26  stations  is  -4  cm.  The  Sea  Surface  Topography 
variation  and  the  location  of  the  tide  gauges  is  shown  on  Figure  4. 
Therefore,  the  value  of  Ah(<f>,  a)  over  the  North  American  Continent 
was  chosen  as  -4  cm.  The  choice  of  the  constant  value  of  the  distortion 
for  the  case  of  the  American  continent  should  be  understood  as  merely 
an  approximation.  Due  to  irregular  distribution  of  SST  and  improper 
numerical  procedure  in  the  1929  adjustment  the  distortion  function 
is  not  constant  but  actually  reflects  an  irregular  internal  warp  in 
the  continental  vertical  datum.  This  will  be  studied  in  the  sequel. 

b)  South  America 

No  information  is  available  on  vertical  datum  or  datums  on  this  continent. 
However,  since  we  are  interested  in  the  general  effect  of  discrepancies 
of  the  order  of  variations  of  the  SST,  it  is  sufficient  to  numerically 
simulate  a  particular  value  of  Ah  in  this  area.  Let  us  assume  that 
the  unified  continental  levelling  net  was  constructed  and  the  adjustment 
was  done  with  a  single  zero-  constraint  at  a  given  tidal  station. 

Suppose  this  tide  gauge  was  located  in  Buenos  Aires  where  the  Sea  Sur¬ 
face  Topography  with  respect  to  280  cm  reference  is  -40  cm  (see  Figure 
10).  In  such  case  the  value  of  Ah(4>,  A)  over  the  S.  American  Con¬ 
tinent  will  be  -40  cm.  The  choice  of  the  master  tide  gauge  in  Buenos 
Aires  is  quite  arbitrary.  From  Figure  10,  we  see  that  if  we  chose 
any  other  location  along  the  coastal  line  we  would  get  the  values  ranging 
from  -80  cm  to  +20  cm. 

c)  Europe 

The  adjustment  of  the  United  European  Levelling  Net  was  carried  out 
with  respect  to  a  single  fixed  datum  point  in  the  Netherlands:  Normal 
Amsterdam  Pei 1  (NAP).  At  this  tide  station  the  sea  Surface  Topography 
with  respect  to  280  cm  reference  is  -55  cm  (Figure  8).  Therefore, 
the  value  of  Ah($,  a)  over  Europe  could  be  taken  as  -55  cm  (Lisitzin, 
1974,  p.  163). 

d)  Africa 

There  is  no  information  on  the  levelling  net  adjustment  on  this  con¬ 
tinent.  Therefore,  assume  the  unified  continental  levelling  net  was 
adjusted  with  a  single  zero  constraint  at  the  tide  gauge  at  Lagos. 

The  Sea  Surface  Topography  at  Lagos  with  respect  to  280  cm  reference 
is  -10  cm  (Figure  11).  Therefore,  the  Ah(<t>,  a)  can  be  taken  as  -10 
cm  over  the  whole  African  Continent.  The  choice  of  master  tide  gauge 
at  Lagos  is  again  arbitrary.  If  we  pick  at  random  a  point  on  the  coast 
line  we  would  get  the  values  for  the  distortion  function  ranging  from 
-60  cm  to  about  +60  cm  (Figure  9). 
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e)  Asia 

There  is  no  information  on  the  levelling  net  adjustment  on  this  continent. 
There,  assume  the  unified  continental  levelling  net  was  adjusted  with 
a  single  zero-constraint  at  Leningrad.  The  Sea  Surface  Topography 
at  Leningrad  with  respect  to  280  cm  reference  is  approximately  -30  cm 
(Figure  9).  Therefore,  the  value  of  a h (<J> »  x)  over  Asia  will  be  taken 
as  -30  cm.  If  we  pick  any  other  location  for  the  master  tide  gauge 
we  would  obtain  different  values  for  the  distortion  function  ranging 
from  about  -120  cm  to  about  +80  cm. 

f)  Australia 

The  adjustment  of  the  Australian  Height  Datum  (1971)  was  carried  out  with 
the  heights  fixed  to  zero  at  30  tide  gauges  around  the  coast  of  the 
continent.  The  Sea  Surface  Topography  average  over  those  tide  gauges 
(with  respect  to  280  cm  reference)  is  +35  cm  (Figure  6).  Therefore, 
the  value  of  AhU,  x)  over  Australia  could  be  taken  as  +35  cm. 

On  the  other  hand,  the  free  network  has  been  calculated  (Mitchell, 
1972)  with  a  single  constraint  at  Jervis  Bay,  where  the  SST  is  equal 
to  30  cm.  Since  a  single  constraint  approach  is  more  accurate  than 
the  overconstrained  problem  (as  was  shown  in  Chapter  4)  we  will  use 
the  value  of  +30  cm  to  represent  Ah  at  this  area. 

g)  Antarctica 

There  is  no  information  on  the  vertical  datum  on  this  continent,  neither 
is  there  on  Sea  Surface  Topography  below  the  latitude  60  South.  There¬ 
fore,  Antarctica  will  be  treated  as  the  ocean  area.  The  constant  value 
of  Ah(<f>,  x)  over  Antarctica  will  be  taken  as  zero  cm. 

In  order  to  study  separately  the  effects  of  different  components 
of  the  height-error  function  Ah(4>,  x),  we  constructed  three  different 
global  models  of  this  function  a!^  ,  Ah2,  Ah3,  proceeding  from  the 
simplest  to  the  more  complicated  one. 

Numerically,  all  our  models  are  constructed  as  a  set  of  64800 
mean  values  based  on  the  l°x  1°  blocks  covering  the  earth's  surface 
in  a  regular  equiangular  grid.  The  continental  outline  with  the 
l°x  1°  resolution  was  obtained  based  on  the  set  of  l°x  1°  mean  land 
elevations  and  ocean  depths  provided  by  the  Defense  Mapping  Agency 
Aerospace  Center  in  1979  (archived  as  file  no.  11  on  the  GS160  mag¬ 
netic  tape.  Department  of  Geodetic  Science  and  Surveying,  The  Ohio 
State  University). 

Figure  17  shows  the  three-dimensional  view  of  function  Ahj  which 
is  a  simple  patchwise  step  function  assuming  constant  value  over  different 
continental  plates: 
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The  reasoning  for  such  values  was  given  above.  In  the  ocean  areas 
our  step  function  Ahi  is  set  to  zero,  which  amounts  to  the  assump¬ 
tion  that  the  mean  sea  level  coincides  with  our  model  geoid.  The  pur¬ 
pose  of  this  first  model  is  to  reflect  the  relative  displacements  between 
different  continental  vertical  datums.  The  individual  vertical  datums 
are  considered  here  the  undisturbed  flat  reference  surfaces.  The  range 
of  Ahi  values  used  in  this  numerical  model  is  from  -55  to  +30  cm. 

In  our  second  model-function  Ah2  ,  we  introduce  the  SST  variation 
over  the  ocean  regions.  The  continental  values  of  Ah2  are  exactly 
the  same  as  in  the  first  model.  The  purpose  of  this  model  is  to  reflect 
the  relative  displacements  of  different  vertical  datums,  understood 
as  the  rigid  vertical  shiftings  between  the  undisturbed  continental 
reference  planes,  against  the  slowly  varying  SST  over  the  oceans. 

Figure  18  shows  the  3-dimensional  view  of  the  model  2  error  function. 

The  magnitude  of  Ah2  varies  within  the  interval  (-160  cm,  +120  cm). 

The  values  of  the  Sea  Surface  Topography  were  generated  using  the  sub¬ 
routine  LISITZ  by  Kostas  Katsambalos,  of  the  Ohio  State  University 
in  1977.  This  subroutine  uses  the  digitized  values  of  SST  from  the 
Lisitzin’s  world  chart  (Lisitzin,  1974,  Figure  1)  on  the  10°x  10°  grid. 
The  method  of  linear  interpolation  is  used,  and  the  Lisitzin's  global 
average  of  280  dyn.  cm  was  removed.  Katsambalos  (ibid.)  estimates 
the  accuracy  of  the  routine  as  +10  dyn. cm.  The  subroutine  LISITZ  gives 
values  of  SST  only  between  Latitudes  +60  and  -60.  Therefore,  beyond 
this  range  we  assume  the  zero  value  to  represent  our  Ah2  function 
(Figure  18.). 

The  third  and  most  complex  model  is  shown  in  Figure  19.  It  differs 
from  model  2  function  only  in  the  values  over  the  N.  Americal  continent. 
Instead  of  the  constant  reference  level,  the  vertical  datum  distortion 
function  d  introduced  in  Section  4.3  was  used  to  represent  Ah3  over 
this  continent.  The  values  of  the  distortion  function  d  were  con¬ 
structed  based  on  the  grid  of  15x15=225  values  shown  on  Figure  15. 

Between  grid  points,  the  l°x  1°  values  were  estimated  using  the  six- 
point  weighted  moving-average  prediction  method  (Davis,  1973).  The 
purpose  of  this  third  model  is  to  reflect  simultaneously  all  three 
possibile  sources  of  vertical  datum  inconsistencies:  the  relative 
rigid  displacements  between  the  continental  flat  reference  levels, 
the  SST  effect  over  oceans,  and  the  internal  vertical  datum  distortion 
due  to  improper  numerical  procedure  during  the  levelling-data  adjustment. 
The  range  of  Ah3  values  falls  in  the  interval  (-160  cm,  +120  cm). 

We  can  expect,  that  all  three  components  of  the  height  inconsistency 
function  will  contribute  to  the  errors  in  various  gravimetric  quantities. 
The  next  section  will  discuss  that  problem. 


.  THE  EFFECT  OF  THE  INCONSISTENCIES  IN  THE  GLOBAL  VERTICAL  DATUM 


ON  THE  DETERMINATION  OF  SPHERICAL  HARMONIC  COEFFICIENTS  OF  THE  GEOPOTENTIAL 


FROM  GRAVITY  DATA  -  THE  MODELLING  STUDIES 


In  this  section  we  evaluate  the  effects  of  vertical  datum  inconsisten¬ 
cies  on  spherical  harmonics  coefficients  of  the  Earth's  gravitational 
potential.  The  three  different  models  of  height  inconsistency  function 
(developed  in  the  previous  chapter)  will  be  used  as  the  basis  for  this 
study. 

6.1  The  Spherical  Harmonic  Analysis  of  Different  Models  of  the  Height 


Inconsistency  Function  Ah 

Pi rst,  we  compute  the  set  of  corrections  to  the  spherical  harmonic  coef 
ficients  of  Earth's  gravitational  field  based  on  the  spectral  analysis 
of  the  three  model  functions  Ahlt  Ah2,  and  Ah3. 

Consider  the  spherical  harmonic  expansion  of  Earth's  gravitational 
field  in  the  form 

v  =  ~r  (1  +  l  (f)0  l  (C  cos(mx)  +  S  sin(mx))P (sin  *))  (7) 


where  (Heiskanen  and  Moritz,  1967): 


JJ  Ag*  Pnm(sin  $) 

a 


cos(mx) 


sin(mx) 


average  value  of  gravity  (979.8  gals;  Gravity  Formula  1980) 
mean  gravity  anomaly  in  blocks  whose  size  is  da. 
geocentric  gravitational  constant 

the  distance  from  the  coordinate  origin  to  the  point  of 
evaluation 

a  nominal  earth  radius 


longitude  and  geocentric  latitude  of  the  point  of  evaluation 

fully  normalized  associated  Legendre  function  of  degree  n 
and  order  m. 


In  formula  (7)  we  assume  that  the  coordinates  origin  coincides  with 
the  center  of  mass  of  the  earth  (no  first  degree  term  in  the  expansion). 

Values  of  the  potential  coefficients  are  most  currently  found 
from  the  analysis  of  satellite  data  or  through  a  combination  of  gravi¬ 
metric  and  satellite  data  (Rapp,  1973).  However,  we  assume  here  that 
the  coefficients  are  computed  from  terrestrial  gravity  anomalies  using 
formula  (8). 

Let's  assume  Ag*  refers  to  the  idealized  geoid  -  for  the  purpose 
of  this  modelling  study  it  will  be  taken  as  one  implied  by  the  Lisitzin's 
average  sea  surface  topography  of  280  dyn.cm.  Equation  (8)  produces 
coefficients  centered  with  respect  to  a  reference  field  implied  by 
the  normal  potential  of  the  level  ellipsoid  which  most  closely  approxi¬ 
mates  our  idealized  geoid. 

Actually,  instead  of  the  theoretical  Ag*  ,  we  have  access  only 
to  anomalies  Ag  which  were  reduced  from  observed  surface  anomalies 
to  a  'geoid'  in  a  best  possible  way.  The  problem  is,  that  such  a'geoid' 
is  not  unique  but  is  implied  by  the  zero-level  of  a  local  vertical 
datum.  Since  the  zero-height  is  established  in  connection  with  mean 
sea  level  surrounding  a  given  continent,  and  since  mean  sea  level  differs 
from  the  equi potential  surface  of  the  true  geoid  by  irregular  Sea  Surface 
Topography,  different  vertical  datums  in  the  real  world  define  different 
equipotential  surfaces  of  the  actual  gravitational  field  V  .  In  other 
words,  the  zero-height  implied  by  different  vertical  datums  is  not  unique 
but  varies  from  datum  to  datum  by  the  amount  probably  in  the  order  of 
sea  surface  topography  variation.  From  the  above,  we  can  conclude  that 
available  global  mean  gravity  anomaly  data  Ag  should  be  further  reduced 
to  our  model  goid  by  the  height  difference  between  the  zero-height  im¬ 
plied  by  each  vertical  datum  and  the  zero-height  implied  by  our  model 
geoid.  Let  us  denote  this  height  difference  by  Ah  .  Since  for  most  of 
the  world  too  little  information  is  available  on  actual  vertical  datums, 
their  interrelation  with  respect  to  our  model  geoid  can  not  be  precisely 
determined.  Therefore,  for  the  purpose  of  this  study  we  will  use  the 
three  models  of  our  discrepancy  function  constructed  in  the  previous 
chapter.  Functions  Aht,  Ah2,  Ah3  will  serve  as  the  three  different 
approximations  to  the  actual  height-inconsistency  function  Ah  . 

Therefore,  the  results  obtained  should  be  considered  as  an  approximation 
to  the  actual  impact  of  the  true  vertical  datum-error  on  geopotential. 
Suppose  we  had  been  given  gravity  anomalies  Ag  which  are  considered 
to  be  reduced  to  the  geoid  implied  by  an  individual  vertical  datum  D  . 
Originally,  Ag  were  computed  from  the  formula  (Heiskanen,  Moritz,  1967): 


=  90bs 


-  Ifl  h  - 
ah  nD 


where 


observed  value  of  gravity 
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is.  . 

ah 


vertical  gradient  of  gravity 


hp  -  height  difference  between  the  surface  point  of  the  measurement 
and  corresponding  point  on  the  geoid  (implied  by  the  vertical 
datum  D). 

y  -  normal  gravity  on  the  reference  ellipsoid 

On  the  other  hand,  introducing  the  height  correction  Ah  from  the 
particular  local  geoid  implied  by  datum  D  to  our  model  geoid  implied 
by  the  Lisitzin's  SST  global  average,  we  are  interested  in  the  more 
accurate  reduction  of  gravity  data  given  by  tne  formula: 

Ag*  =  9obs  '  fh  (hD  *  Ah)  ”  Y  (9) 


where  Ah  is  the  height  inconsistency  function.  The  Ag*  values 
obtained  from  (9)  should  be  used  in  eq.  (8). 

Now  we  can  compute  the  change  in  gravity  anomaly  caused  by  the 
introduction  of  our  additional  height  correction  Ah: 

6Ag  =  Ag*  -  Ag  =  -yj|  Ah  =  0.3086  Ah  (mgal)(where  Ah  in  meters) 

which  has  the  form  of  the  free-air  correction.  Here  Ah  represents 
the  values  of  our  height  correction  function  Ah(<f%  x),  and  could  be 
approximated  by  our  model-functions  Ahj ,  Ah2,  Ah3. 

Now,  we  can  compute  the  corrections  6Cnm,  ^nm  to  potential 
coefficients  implied  by  the  assumed  inconsistencies  in  vertical  datums 
Ah: 


/  _  \ 

i  sCnm 

cos  mx 

s  °l‘?/86E n5  JJ  ah(*.  a  )  P  (sin  +) 

|  4irG(n-l)  o  nm 

sin  mx 

(  ^nm 

where  Ah(<f>,  x)  (in  meters)  is  height  inconsistency  function  described 
in  Section  5.  The  other  symbols  are  explained  after  formula  (8). 

Instead  of  the  theoretical  function  Ah  the  approximations  Ah^ 
Ah2,  and  Ah3  described  in  Chapter  5  were  used.  The  actual  integration 
was  carried  out  using  Fast  Fourier  Transform  method  implied  in  the  FORTRAN 
subroutine  HARMIN  described  by  Colombo  (1981).  The  subroutine  was 
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run  using  integrated  associated  Legendre  functions  for  geocentric  lati¬ 
tudes,  and  the  optimum  smoothing  factors  recommended  by  Colombo.  The 
advantage  of  HARMIN  is  its  very  efficient  calculation  procedure.  Using 
this  program  a  discrete  spectral  analysis  up  to  degree  180  was  performed 
on  the  three  different  sets  of  64800  l°x  1°  block-values  representing 
our  model  height-inconsistency  functions  Ahi ,  Ah2,  and  A_h 3 .  As  a 
result  three  different  sets  of  correction  coefficients  acnm>  aSnm 
were  generated. 

The  magnitude  of  the  computed  correction  coefficients  was  tested 
against  the  magnitude  of  potential  coefficient  obtained  by  Rapp  (1981) 
from  the  combination  of  the  terrestrial  gravity  data,  SEASAT  altimeter 
data,  and  other  satellite-derived  data  (see  formula  (8)).  The  degree 
variances  of  correction  coefficients  up  to  n=180  were  evaluated  according 
to  the  formula: 


I  (6Cnm  +  5SL> 

nm  nm 


The  relative  magnitude  of  corrections  with  respect  to  the  reference 
coefficient  (Rapp,  1981)  can  be  expressed  by  the  ratio: 


where. 


=  l  (C2  +  S2  ) 

nm  nm' 
m-U 


represents  the  degree  variances  of  reference  coefficients  C  ,  S„  . 
r  3  nm  nm 

Figures  20,  21,  and  22  show  the  relative  magnitude  of  correction 
coefficients  Rn  obtained  for  our  model  datum-inconsistency  functions 
Ahx ,  Ah2,  and  Ah3  respectively.  In  all'  three  cases  the  correction 
coefficients  bear  significant  contribution  to  the  original  coefficients 
only  up  to  degree  60.  Above  degree  60  the  spectrum  flattens,  and  the 
individual  relative  errors  do  not  exceed  the  magnitude  0.002  for  the 
case  of  Ah2  and  Ah3  ,  and  0.001  for  Ah^  We  conclude  that  the 
effect  of  vertical  datum  inconsistency  is  mostly  of  long  wavelength 
nature,  and  that  for  higher  frequencies  the  relative  error  in  the  co¬ 
efficients  is  on  the  order  of  the  background  noise.  The  distribution 
of  the  relative  errors  in  the  frequency  domain  for  models  Ah2  and 
Ah3,  as  shown  on  Figures  21  and  22  respectively,  is  very  much  alike. 


*•  %  «V*  .s  .*•.  *  , 


•  *.’•  w'*  /•  .*• 
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Figure  20.  Relative  magnitude  of  error  coefficients  (by  ratio  of 
degree  variances)  computed  for  the  vertical  datum  in¬ 
consistency  model  Ahi  (see  eq.  12). 


although  the  inclusion  of  the  internal  distortion  of  the  North  American 
datum  (function  Ah3)  caused  slightly  higher  relative  errors  for  most 
degree  variances  (see  also  Table  1). 


Table  1  shows  the  comparison  between  the  magnitude  of  relative 
errors  in  potential  coefficients  expressed  by  formula  (12)  for  three 
different  models  of  vertical  datum  inconsistencies  represented  by 
functions  Ah3,  Ah2,  and  Ah3.  The  maximum  individual  relative  error 
(by  degree)  in  coefficients  for  the  simplest  model  represented  by 
Ahi  (no  SST  included)  occurs  at  degree  4  and  has  the  magnitude  of 
0.00253  (0.253%).  For  the  other  two  models  which  include  SST  the  max¬ 
imum  relative  error  is  about  0.009  (0.9%)  and  occurs  at  degree  6. 
Comparing  the  effects  of  Ah3  (simple  step-function,  no  SST)  with  those 
of  Ah2  or  Ah3  (SST  included)  we  conclude  that  the  variation  of 
SST  over  oceans  brought  about  a  significant  increase  of  power  spectrum 
of  correction  coefficients.  The  models  containing  SST  produce  the 
spherical  harmonics  coefficients  that  are  much  larger  in  magnitude 
than  those  produced  by  a  simple  step-function  Ahi  .  The  addition 
of  Sea  Surface  Topography  had  a  very  large  effect  of  magnifying  the 
errors  more  than  3  times.  The  oceans  cover  more  than  70%  of  our  globe, 
and  therefore  the  variations  over  oceans  determine  the  overall  character 
of  the  inconsistency  model.  If  we  multiply  (12)  by  100  we  obtain  the 
average  percentage  correction  to  potential  coefficients  due  to  the 
datum  inconsistency.  From  Table  1  we  conclude  that  the  global  datum 
inconsistency  can  produce  a  significant  error  in  the  determination 
of  potential  coefficients.  For  low  degree  harmonics  the  relative  error 
can  be  on  the  order  of  1  percent.  By  significant  we  understand  the 
error  of  the  relative  magnitude  reaching  1%  level. 

The  spectral  analysis  of  our  three  models  of  the  spatial  incon¬ 
sistencies  in  the  global  vertical  datum  (as  represented  on  Figures 
20,  21  and  22)  shows  significant  effect  only  on  the  low  degree  harmonics. 
This  is  a  result  of  the  specific  nature  of  our  approximating  functions 
Ahi,  Ah2,  and  Ah3  (see  Figures  17,  18  and  19).  All  of  them  represent 
a  flat  or  slowly  varying  surface,  showing  no  local  features  of  short- 
wavelength  characteristics  (except  perhaps  discontinuities  along  con¬ 
tinental  edges).  Therefore,  the  particular  distribution  of  energy 
in  the  spectrum  could  have  been  expected.  We  should  mention  also  that 
the  digital  model  of  SST  used  in  the  computations  was  based  on  the  rather 
crude  10°x  10°  grid,  and  this  certainly  could  add  up  to  the  overall 
smoothing  of  the  distortion  spectrum  (high-cut  filter  effect). 

More  than  two  thirds  of  the  entire  energy  in  the  spectrum  comes 
from  the  discrepancy  between  MSL  and  the  geoid  over  the  oceans,  as 
could  be  seen  by  comparing  the  first  column  with  the  second  or  third 
columns  in  Table  1.  The  relative  error  due  to  the  internal  distortion 
in  the  N.  American  vertical  datum  contributes  very  little  to  the 
cumulative  effect;  however,  it  may  be  significant  in  some  applications. 

6.1.1  The  Error  (by  Degree  Variances)  in  Geoid  Undulations,  Gravity 
Anomalies,  and  Deflections  of  Vertical  due  to  Vertical  Datum  Incon- 
sistency  Models 
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Table  1 

The  Relative 

Error  in  Potential 

Coefficients  (by  Ratio  of  Square 

Root  Degree  Variances,  eq. 

12)  due  to  the  three  Vertical 

pH 

Datum  Inconsistency  Models:  Ah,,  Ah2, 

and  Ah3. 

(Unitless  quantities) 
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‘.\v 
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DEGREE (N) 
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Ah2 

Ah3 
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2 

3.62E-05 

1.37E-04 
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3 

2.11E-03 

4.76E-03 
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8 

1.79E-03 
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5.126-03 

9 
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4.316-03 

4.426-03 
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10 

1.83E-03 

5.22E-03 

5.26E-03 

11 

1. 72E-03 
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12 
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7.75E-03 
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13 

1.42E-03 

3.86E-03 
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i>  A 

14 

1.90E-03 

4.82E-03 

4.92E-03 

oil 

15 

2.376-03 

5.75E-03 
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16 

1.43E-03 

5.70E-03 
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17 

2.00E-03 

5.30E-03 

5.236-03 
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1.61E-03 
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5. 146-03 
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19 

2.00E-03 

4.58E-03 

4.566-03 

20 

2.35E-03 

4.19E-03 

4.195-03 

21 

2.05E-03 

3.81E-03 
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22 

1.43E-03 

4.98E-03 

5.02E-03 
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23 

1.796-03 
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i  *  •  • 
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where  R  denotes  the  mean  earth  radius.  The  contribution  to  the  geoid 
undulation  in  terms  of  the  square  roots  of  undulation  degree  variances 
due  to  our  three  models  of  the  inconsistencies  in  global  vertical  datum 
is  shown  in  Figures  23,  24  and  25  respectively.  In  Table  2  we  compare 
the  effects  (up  to  degree  50)  on  geoid  undulation  produced  by  our  three 
models  of  vertical  datum  inconsistencies  Ah^  Ah2  and  Ah3.  The 
maximum  square  root  of  undulation  degree  variance  was  computed  as  11.16  cm 
(at  degree  2)  for  the  first  model,  42.38  cm  (at  degree  2)  for  the  second 
model  and  42.05  cm  (at  degree  2)  for  our  third  model. 


The  existence  of  the  zero  degree  unulation  term  in  the  error  spectrum 
can  be  interpreted  in  the  following  way.  Theoretically,  if  we  evaluated 
the  anomalous  gravity  potential  T  with  respect  to  the  reference  normal 
potential  of  the  level  ellipsoid,  which  is  chosen  to  have  a  mass  equal 
to  the  mass  of  the  Earth  and  potential  equal  to  that  associated  with 
the  geoid,  then  the  zero-degree  coefficient  in  the  spherical  harmonic 
expansion  of  T  or  undulation  N  would  vanish.  This  is  true  only 
if  the  terrestrial  gravity  data  used  were  properly  reduced  to  the  geoid. 
The  geoid  is  approximated  by  a  given  level  ellipsoid  which  induces 
normal  gravity  used  in  the  gravity  data  reduction.  Geometrically  this 
approximation  means  that  the  average  (over  the  unit  sphere)  deviation 
of  the  geoid  from  the  reference  ellipsoid  is  zero,  and  the  origin  of 
the  ellipsoid  coincides  with  the  center  of  mass  of  the  Earth. 

In  practice  however  we  deal  with  the  gravity  data  that  have  been 
reduced  not  to  the  idealized  geoid,  but  to  some  nearby  surface  which 
is  not  necessarily  the  equipotential  surface.  Let  us  for  the  moment 
call  this  surface  a  pseudogeoid.  The  height  error  function  Ah 
described  in  the  previous  chapter  represents  precisely  the  deviation 
of  pseudogeoid  from  the  geoid.  In  the  frequency  domain,  this  de¬ 
viation  is  represented  by  the  error  power  spectrum  xn,  that  is  by 
the  various  harmonic  components  of  the  difference: 


{pseudogeoid  undulation}  -  {geoid  undulation}  :=  6N 


The  non-vanishing  zero  degree  undulation  term  in  the  error  spectrum 
shows  that  the  average  height  of  the  pseudogeoid  with  respect  to  our 
reference  normal  ellipsoid  does  not  reduce  to  zero.  The  detailed  inter¬ 


pretation  of  the  zero-degree  term  can  be  found  for  example  in  Heiskanen 
and  Moritz  (1967).  In  geometric  terms  the  semimajor  axis  of  the  mean 
best  fit  ellipsoid  associated  with  the  pseudogeoid  differs  from  the 
semimajor  axis  of  the  normal  ellipsoid  associated  with  the  geoid  by 
the  amount  expressed  by  the  zero  degree  term  in  the  error  spectrum. 

By  the  best-fit  ellipsoid  implied  by  the  pseudogeoid  we  understand  the 
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Table  2 

The  Error  in  Geoid  Undulation  (by  Square  Root  Degree  Variances  up  to 
Degree  50,  see  eq.  13)  due  to  the  three  Vertical  Datum 
Inconsistency  Models:  Ahj ,  Ah2,  and  Ah3. 

Units  are  cm. 


REE(N) 

Ah  j 

CM 

.C 

<1 

Ah3 

0 

10.562 

3.369 

3.738 

2 

11. 162 

42.384 

42.051 

3 

3.996 

9.008 

9.239 

4 

2.562 

8.241 

8.229 

5 

1.350 

4.406 

4.453 

6 

1.218 

5.192 

5.204 

7 

0.675 

1.957 

1.970 

8 

0.551 

1.554 

1.580 

9 

0.582 

1.  167 

1.  195 

10 

0.414 

1.  179 

1.188 

11 

0.292 

0.969 

0.962 

12 

0.258 

0.818 

0.828 

13 

0.206 

0.560 

0.551 

14 

0.164 

0.416 

0.425 

15 

0.  191 

0.464 

0.471 

16 

0.136 

0.540 

0.541 

17 

0.  151 

0.401 

0.399 

18 

0.112 

0.356 

0.358 

19 

0.125 

0.286 

0.285 

20 

0.111 

0.198 

0.198 

21 

0.104 

0.  194 

0.196 

22 

0.086 

0.298 

0.300 

23 

0.082 

0.217 

0.216 

24 

0.064 

0.193 

0.194 

25 

0.068 

0.187 

0.186 

26 

0.065 

0.119 

0.123 

27 

0.058 

0. 139 

0.140 

28 

0.063 

0.  184 

0.184 

29 

0.055 

0.  132 

0.133 

30 

0.062 

0.  163 

0.163 

31 

0.048 

0.113 

0.114 

32 

0.04  9 

0.090 

0.090 

33 

0.048 

0.084 

0.083 

34 

0.047 

0.131 

0.131 

35 

0.042 

0.109 

0.110 

36 

0.041 

0.098 

0.098 

37 

0.040 

0.080 

0.080 

38 

0.040 

0.057 

0.056 

39 

0.041 

0.064 

0.065 

40 

0.029 

0.075 

0.076 

41 

0.031 

0.070 

0.070 

42 

0.030 

0.081 

0.081 

43 

0.027 

0.064 

0.064 

44 

0.030 

0.054 

0.054 

45 

0.025 

0.047 

0.047 

46 

0.024 

0.064 

0.064 

47 

0.024 

0.055 

0.055 

48 

0.025 

0.065 

0.066 

49 

0.023 

0.055 

0.055 

50 

0.019 

0.039 

0.039 
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The  error  in  geoid  undulation  (by  square  root  degree 
variances,  eq.  13)  due  to  the  vertical  datum  incon¬ 
sistency  model  Ahj. 

Units  are  cm. 
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The  error  in  geoid  undulation  (by  square  root  degree 
variances,  3q.  13)  due  to  the  vertical  datum  incon¬ 
sistency  model  Ah 3. 

Units  are  cm. 


Table  3 

The  Error  in  Gravity  Anomalies  (by  Square  Root  Degree  Variances,  eq.  15) 
due  to  the  Three  Vertical  Datum  Inconsistency  ModelsiAhi,  Ah2  and  Ah3. 

Units  are  micro-gals. 
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DEGREE(N)  Aht  Ah2  Ah3 


0 

16.24 

5.18 

5.75 

1 

15.96 

74.57 

74.34 

2 

17.16 

65.18 

64.6  7 

3 

12.29 

27.71 

23.42 

4 

11.62 

38.02 

37.97 

5 

8.3C 

27.  10 

27.39 

6 

9.36 

39.92 

40.02 

7 

6.22 

18.06 

13.13 

e 

5.93 

16.72 

17.00 

9 

7.17 

14.36 

14.70 

10 

5.72 

16.31 

16.45 

11 

4.49 

14.91 

14.80 

12 

4.37 

13.84 

14.00 

13 

3.80 

10.33 

10.16 

14 

3.28 

8.32 

8.49 

15 

4.11 

10.00 

10.14 

16 

3.  14 

12.46 

12.47 

17 

3.72 

9. Bo 

9.83 

18 

2.93 

9.31 

9.35 

19 

3.46 

7.92 

7.89 

20 

3.24 

5.76 

5.79 

21 

3.21 

5.97 

6.01 

22 

2.77 

9.62 

9.70 

23 

2.  77 

7.33 

7.31 

24 

2.25 

6.32 

6.37 

25 

2.51 

6.91 

6.88 

26 

2.49 

4.53 

4.72 

27 

2.33 

5.56 

5.59 

28 

2.60 

7.  of 

7.63 

29 

2.35 

5.69 

5.74 

30 

2.77 

7.29 

7.25 

31 

2.23 

5.23 

5.25 

32 

2.33 

4.28 

4 . 2  o 

33 

2.34 

4.12 

4.10 

34 

2.38 

6.e>3 

6.66 

35 

2.21 

5.71 

5.74 

36 

2.22 

5.27 

5.29 

37 

2.22 

4.43 

4.45 

38 

2.29 

3.22 

3.21 

39 

2.39 

3.  76 

3.77 

40 

1.72 

4.52 

4.58 

41 

1.89 

4.29 

4.29 

42 

1.92 

5.09 

5.11 

43 

1.77 

4. 16 

4.15 

44 

1.95 

3.56 

3.53 

45 

1.69 

3.21 

3.20 

46 

1.67 

4.43 

4.42 

47 

1.70 

3.89 

3.90 

48 

1.80 

4.72 

4.  74 

49 

1.66 

4.03 

4.07 

50 

1.44 

2.91 

2.92 
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Table  4 


The  Error  in  Deflections  of  the  Vertical  (by  Square  Root  Degree 
Variances,  eq.  17)  due  to  the  Vertical  Datum  Inconsistency  Model  Ah3. 

Units  are  Seconds  of  Arc. 


DEGRtE(N) 


3 

4 

5 

6 
7 
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10 

19 
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21 
22 
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24 
23 
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34 
33 
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30 

39 

40 

41 

42 

43 

44 
43 

46 

47 
40 

49 
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3.332-02 
i . C4E-02 
1 . 192-02 
7. C0E-03 
1 . 092-C2 
4 . 772-03 
4 . 342-03 
3. 672-03 
4.032-C3 

o  rr 
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ellipsoid  producing  no  zero-degree  undulation  term  (Heiskanen,  Moritz, 
1967,  eg.  (2-194a)).  In  geometric  terms,  the  pseudogeoid  encloses 
a  slightly  different  volume  than  the  geoid  does.  In  other  words  the 
volume  of  a  thin  irregular  layer  described  by  our  height  error  function 
&h  does  not  average  to  0. 

This  does  not  necessarily  contradict  the  fact  that  the  average 
SST  used  as  a  primary  component  in  Ah  was  claimed  to  be  zero. 

First  of  all  the  Lisitzin's  280  dyn.cm  average  of  SST  variation  was 
implied  by  data  irregularly  distributed  between  latitudes  +60°  and  -60°. 
Near  the  poles  there  was  no  data  at  all. 

Secondly,  the  280  dyn.cm  average  pertains  only  to  the  oceanic 
part  of  Earth's  surface,  whereas  our  model  height  error  function  covers 
the  entire  globe. 

In  the  next  step  of  this  analysis,  for  each  model  we  computed  the 
cumulative  effect  on  geoid  undulation  up  to  degree  180: 


(14) 


This  effect  is  16.24  cm  for  function  Ahj ,  44.90  cm  for  Ah2  and 
44.67  cm  for  model -function  Ah3.  Again  almost  two-thirds  of  the 
total  effect  comes  from  the  SST  variation  over  oceans,  while  the  con¬ 
tribution  from  the  relative  rigid  displacements  of  the  flat  continental 
datums  (function  Ahx)  amounts  only  to  about  16  cm.  The  inclusion 
of  the  internal  distortion  of  the  N.  American  vertical  datum  had  a 
cumulative  effect  on  the  order  of  millimeters. 

The  computed  effect  of  inconsistencies  in  the  global  vertical 
datum  can  also  be  expressed  in  terms  of  the  gravity  anomaly  degree 
variances: 


where  G  :=  the  average  value  of  gravity. 
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Figure  27.  The  error  in  gravity  anomaly  (by  square  root  degree  variances, 
qq.  15)  due  to  vertical  datum  inconsistency  model  Ah2. 

Units  are  mgals. 
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Figure  28.  The  error  in  gravity  anomaly  ( by  square  root  degree  variances, 
eq.  15)  due  to  the  vertical  datum  inconsistency  model  Ah3. 
Units  are  mgals. 
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Figures  26,  27,  and  28  show  the  contribution  to  gravity  anomaly 
by  the  square  root  of  cn  (formula  (15))  produced  by  model  1,  2,  and 
3.  Table  3  compares  the  errors  in  the  gravity  anomaly  represented 
by  the  above  three  figures  for  our  model  functions  Ahj ,  Ah2  and  Ah3 
respectively.  The  maximum  individual  errors  (by  degree)  reach  the 
magnitude  0.017  mgal  at  degree  2  for  the  case  of  the  simplest  model 
Ah3  ,  up  to  0.074  mgal  (at  degree  1)  for  the  most  complex  one.  Also, 
we  computed  the  cumulative  error  in  gravity  anomaly  expressed  by  the 
formula: 


180  n 

i  G2("-I)2Jn  +  <16> 

n-0  m=0 


This  error  is  0.043  mgal  for  the  simplest  model  Ahi ,  0.13  mgal  for 
model  Ah2  and  0.13  mgal  for  the  most  complicated  model.  Again  the 
dominant  error  comes  from  the  Sea  Surface  Topography  variation  over 
oceans.  The  effect  of  distortions  in  N.  Americal  vertical  datum  is 
only  on  the  order  of  IE-4  mgal. 

We  also  computed  the  effect  of  the  global  vertical  datum  incon¬ 
sistency  in  terms  of  deflections  of  vertical  degree  variances: 

n 

0  *s  n(n+l)  £  ( <5C^  +  <sS^_)  (17) 

n  '  L  nm  nm 

m=0 


Figure  29  shows  the  frequency  spectrum  decomposition  of  the  error  in 
the  deflection  of  vertical  (by  degree)  generated  by  our  most  complex 
model -function  Ah3,  reflecting  the  rigid  relative  shiftings  of  the 
continental  reference-planes,  SST  over  oceans, and  internal  distortion 
of  the  N.  American  vertical  datum.  The  individual  values  were  computed 
as  square  roots  of  the  deflection  of  vertical  degree  variances  given 
by  formula  (17).  The  units  are  seconds  of  arc.  Also,  in  Table  4  we 
show  those  results  up  to  degree  50.  The  maximum  individual  error  occurs 
at  degree  2  and  has  the  magnitude  of  3.33E-2  sec  of  arc,  whereas 
the  cumulative  effect  up  to  degree  180: 


is  5.2E-2  sec  of  arc. 

i 
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Figure  29.  The  error  in  deflections  of  vertical  (by  square  root  degree 
variances,  eq.  17)  due  to  the  vertical  datum  inconsistency 
model  Ah 3. 

Units  are  seconds  of  arc. 


6.2  Some  Remarks  on  the  Present  Procedures  for  Determination  of  Spherical 
Harmonics  Coefficients  of  Geopotential 

In  the  previous  -section  we  studied  the  effects  of  vertical  datum 
inconsistencies  on  the  determination  of  geopotential,  as  if  the  only 
information  available  was  the  terrestrial  gravity  data.  However,  the 
most  current  procedure  of  computing  the  potential  coefficients  is  through 
a  combination  of  gravimetric  and  satellite-derived  data  (Rapp,  1973). 

Those  procedures  use  satellite  data  to  derive  low  degree  coefficients 
and  terrestrial  data  to  supply  information  for  the  higher  frequency 
part  of  the  spectrum.  Therefore,  only  part  of  the  error  estimated 
in  the  previous  chapter  will  actually  contribute  to  the  error  in  the 
determination  of  the  spherical  harmonic  coefficients.  This  part  comes 
only  from  the  higher  degree  error  coefficients,  produced  by  our  models 
of  inconsistencies  in  the  global  vertical  datum. 

It  might  be  of  interest  to  estimate  the  cumulative  error  in  coef¬ 
ficients  due  to  the  datum  inconsistencies  only  above  certain  degree. 

It  would  reflect  the  fact  that  only  higher  degree  coefficients  are 
affected  by  errors  of  that  nature,  and  that  low  degree  harmonics  are  derived 
in  some  different  way.  However,  the  satellite-derived  coefficients 
are  affected  by  different  sources  of  errors  and  generally  are  less 
accurate  for  higher  degrees.  Therefore,  this  type  of  analysis  could 
provide  an  additional  insight  into  the  nature  of  this  combined  procedure. 
Also,  it  could  suggest  the  specific  degree  beyond  which  the  terrestrial 
information  can  safely  be  incorporated  into  the  combined  solution.  The 
idea  is  to  balance  the  errors  in  the  satellite-derived  coefficients 
with  those  inherent  in  the  terrestrial  data. 

Figure  30  shows  the  cumulative  error  in  the  geoid  undulation  due 
to  our  third  model  (Ah 3 )  of  datum  inconsistency.  The  6  curves  on  this 
figure  represent  the  cumulative  contribution  to  the  undulation  (as 
a  function  of  degree  n)  beyond  the  specified  degree.  This  truncation 
degree  can  be  found  at  the  intersection  with  the  horizontal  axis. 

To  construct  the  curves  we  used  the  formula: 


P  n  i  _ 

CUMULAT.ERROR.lv  (n)  =  R  £  £  (sC?m  +  5S?J  ;  p  <  n  <  180  (18) 

\i=p  m=0 


where  the  starting  degree  p  takes  on  the  values  10,  20,  ...,  60. 

From  Figure  30  we  learn  that  the  cumulative  error  undulation,  when 
the  summation  starts  at  degree  10,  can  amount  to  almost  2  cm  at  degree 
180  (the  highest  curve  on  Figure  30).  If  we  start  the  summation  from 
degree  40  this  cumulative  error  can  reach  0.3  cm  and  for  the  summation 
starting  at  60  it  is  only  0.15  cm  at  degree  180  (the  lowest  curve  on 
Figure  30). 

The  similar  type  of  analysis  for  gravity  anomalies  is  shown  in 
Figure  31.  The  curves  were  constructed  according  to  the  formula: 


Figure  31.  Cumulative  error  contribution  to  gravity  anomalies  (eq.  19) 
beyond  specified  degree  due  tp  the  vertical  datum  inconsis¬ 
tency  model  ah  . 

Units  are  mgal . 
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p  <  n  <  180  for  p  =  10,  20,  ....  60. 


The  correction  coefficients  6Cnm,  6Snm  produced  by  our  third  model 
(Ah3)  of  datum  inconsistency  were  used.  From  the  highest  curve  in 
Figure  31  we  see  that  neqlecting  the  errors  in  the  first  10  degrees 
the  cumulative  effect  (19)  of  errors  at  the  remaining  part  of  the  spectrum 
can  reach  0.053  mgal  at  the  degree  180.  On  the  other  hand  neglecting 
the  errors  of  coefficients  having  degree  less  or  equal  60,  the  cumula¬ 
tive  error  in  the  gravity  anomaly  due  to  datum  inconsistency  beyond 
this  degree  can  only  reach  0.023  mgal  for  n=180  (see  formula  (19)). 

7.  THE  EFFECT  OF  VERTICAL  DATUM  INCONSISTENCIES  ON  THE  DETERMINATION 
OF  GEO ID  UNDULATION  FROM  A  COMBINATION  OF  SPHERICAL  HARMONIC  COEFFICIENTS 


AND  GRAVITY  IN  CAP 

In  practice,  the  geoidal  undulation  are  often  computed  using  Stokes' 
equation  where  the  global  data  set  of  gravity  anomalies  has  been  truncated 
to  form  a  spherical  cap  surrounding  the  point  computation.  Additionally 
this  truncated  terrestrial  data  is  usually  supported  by  the  long  wave¬ 
length  global  information  in  the  form  of  low  degree  harmonic  expansion 
of  geopotential.  The  theory  of  this  combination  of  satellite-derived 
and  terrestrial  data  to  obtain  the  geoid  undulation  (or  disturbing 
potential)  originates  with  M.S.  Molodenskii  (Molodenskii  et  al . ,  1962) 
and  was  further  developed  by  (Heiskanen  and  Moritz,  1967),  (Meissl,  1971), 
(Jekeli,  1980)  and  others. 

The  general  idea  came  from  the  fact  that  global  gravity  data  is 
needed  in  the  classical  Stokes'  equation: 


N  •  4^s  H  Ag  SU)  da 


In  the  equation  the  symbol  S(ij>)  denotes  the  Stokes'  function: 


suj  -  y 


P  (cos  if;) 


i it  -  the  spherical  distance  between  the  variable  surface  element  do 

and  the  point  at  which  the  formula  is  evaluated 

Pn(cos<j>)-  Legendre's  polynomial 

Normally,  if  we  do  not  associate  any  additional  assumptions  with 
Ag  in  equation  (20)  the  left  hand-side  quantity  N  represents  the 
geoid  undulation  from  which  the  zero  and  first  degree  harmonics  have 
been  removed  (Heiskanen  and  Moritz,  1967,  eq.  (2-1636)),  even  if  Ag 
contains  the  zero  and  first  degree  harmonic  components.  In  other  words 
those  two  components  in  Ag  actually  do  not  contribute  to  the  geoid 
undulation  produced  by  eq.  (20). 

Usually,  the  information  on  the  zero-degree  undulation  term  N 
must  be  supplied  to  equation  (20)  in  the  form  of  a  separate  term 


N  =  N°  +  4^G  ^ Ag  da  (21) 


where  (Heiskanen  and  Moritz,  1967,  p.  102): 


N°  =  A9o  +  2GR 


where 

R  -  mean  radius  of  earth 

G  -  mean  value  of  the  normal  gravity  over  the  earth 

<5M  -  mass  difference  between  the  earth  and  reference  ellipsoid 

k  -  gravitational  constant 

Ag0  -  the  zero-degree  component  possibly  present  in  the  spherical 
harmonic  expansion  of  Ag  field. 

If  we  are  using  gravity  anomalies  in  a  spherical  cap  surrounding 
the  point  of  computation,  a  special  correction  term  is  needed  to  account 
for  the  effect  of  truncation  of  the  data  domain.  More  specifically, 
the  gravity  anomaly  data  in  the  cap  are  used  with  the  ordinary 
Stokes'  kernel  SU). 


(23) 
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where 


spherical  cap  around  the  point  of  computation. 


then  the  truncation  error  can  be  expressed  by  (Jekeli,  1980) 


4N‘  "  IS  l,  «in  L% 

n-d 


(24) 


where 


Ag  =  l  Ag 
n=0  r 


denotes  the  expansion  of  the  terrestrial  gravity  anomalies  in  to  the 
series  of  Laplace  spherical  harmonics, 

QlnUc)  =  f  pn(^  sin^  ^ 

c 

(Molodensky  truncation  coefficients),  and 

4>c  -  the  radius  of  spherical  cap. 

If  the  gravity  anomaly  data  in  the  cap  are  used  with  the  modified 
kernel  SU)  -  S(^c): 


4ttG 


II  Ag  (S(cos  })  -  S(cos  4i  ))  da 


(25) 


then  the  truncation  error  can  be  expressed  (Molodenskii  et  al . ,  1962)  as 


6N2  =  H  nl2  Q2n^  Agn 


(26) 
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where 


are  known  as  the  modified  Molodenskii  truncation  coefficients. 

The  Agn  terms  above  are  the  Laplace  surface  harmonics  of  the 
gravity  anomaly  field  Ag.  This  representation  can  be  obtained  from 
a  priori  given  set  of  potential  coefficients  according  to  the  formula: 


n 

Agn  =  G(n_1)  J=Q  {Cmc05  mx  +  Snmsin  mA  }  Pnm(cos  9  }  (27) 


where 

Cnm’  ^nm  “  anomalous  potential  coefficients  relative  to  the  potential 
field  implied  by  the  assumed  reference  ellipsoid 

A*  9  -  polar  coordinates  of  the  point  of  computation 

The  actual  computational  procedure  for  the  evaluation  of  geoid 
undulation  based  on  the  above  method  is  given  by  Rummel  and  Rapp  (1976, 
eq.  21)  (excluding  atmospheric  effect)  as: 


N  =  N°  +  4^G  M  SU)  Ag  do  +  4  l  Q.nUJ 


c  '  2G  '*in"V  A% 


(28) 


where 


No  =  No  "  *Tg  H  S(^  A9o  doc 


(29) 


Substituting  for  N  (eq.  22)  and  using  the  identity  (Jekeli,  1980, 
eq.  123,  124): 


"4?G  H  S{^  A90  doc  =  ^  Qio  Ago 
•  °c 

we  rewrite  equation  (29)  in  more  explicit  form: 


(30) 


/  V 


(31) 


W  —  6^  R  .  /  i  r\  \ 

N°  "  2GR  ‘  2G  Ago  (1'Qio) 

In  order  to  evaluate  the  contribution  of  Ag0  term  to  N0  we  need 
values  of  Q10  and  Q2o  ( Q 2 o  mus t  be  used  in  eq.  (31)  if  the  modified 
kernel  was  used  in  Stokes'  integral). 

It  can  be  shown  that  (Molodenskii  et  al . ,  1962) 


Qi0(*c>  =  _4t  +  5t2  +  6t3  "  7t4  +  (6t2  -  6t4)  lnt( 1+t)  (32) 

where 

.  . 

t  =  sin 


Also,  it  can  be  shown  that  (Jekeli,  1980  eq.  65) 


Q20(l»c)  =  Q10  +  s^c}  (1  '  cos  *c}  (33) 

In  the  previous  section  we  noted  that  the  zero-degree  harmonic 
component  in  the  error  in  gravity  anomaly  field  due  to  the  inconsistency 
in  the  vertical  datum  can  reach  the  magnitude  (Table  3): 


AgQ  =  0.57  x  10-2  mgal  (34) 

This  error  component,  in  turn,  will  be  responsible  for  the  truncation 
effect  in  the  zero-degree  undulation  term  when  the  anomalies  in  the 
cap  are  used  in  the  equation  (28)  with  unmodified  kernel  S(^)  or  with 
modified  one:  SU)  -  S(^c). 

First  we  define  the  truncation  effects  on  zero  degree  undulation 
term  which  contribute  to  the  quantity  defined  by  equation  (31)  as: 

5N?-a«1o4S0  <35> 


or 

6N°  =  £  Q  Ag 
2  2G  y20  yo 

depending  on  the  kernel  used  in  the  Stokes'  integral 

-72- 


(36) 


Figure  32  shows  the  absolute  values  of  these  effects  plotted  as 
functions  of  the  truncation  angle.  We  see  that  for  small  caps  (trun¬ 
cation  radius  ipc  <  38°)  the  modified  Molodenski  kernel  (25)  produces 
smaller  truncation  effect  than  the  unmodified  kernel  (23).  For  the 
cap  radius  ipc  =  38  degrees  both  methods  produce  the  same  truncation 
effect  of  about  2.1  cm.  If  we  consider  larger  caps  (<pc  >  38°)  the 
truncation  effect  of  modified  kernel  is  generally  speaking  greater 
than  that  associated  with  unmodified  one,  and  can  reach  3.6  cm  for 
the  cap  72°  in  radius  and  11.4  cm  for  the  full  cap  180°  in  radius. 

This  last  number  has  no  practical  meaning  since  the  global  set  of  data 
is  required  so  in  fact  there  is  no  truncation  at  all.  It  merely  re¬ 
flects  the  nominal  limit  of  the  expression  (36)  when  ipc  approaches 
180°.  We  rewrite  equation  (36)  in  the  equivalent  form  (Jekeli, 
1980): 


6N°  = 

2 


where 


'4?G  Ag  H  S(c0s  ^  da  +  4^G  Ag  S(^c)  M  do 

a  a 

c  c 

denotes  the  spherical  cap  of  radius  <kc. 


We  see  that  in  the  limit  for  at  -*■  ai80°:=  a  (the  area  of  the 
unit  sphere)  the  integral  in  the  first  term  goes  to  0  from  the  orthogonal 
properties  of  Legendre  polynomials  and  the  integral  in  the  second  term 
goes  to  4rr  which  is  the  area  of  the  unit  sphere.  Therefore,  for 
ipc=180°  we  get  the  nominal  value 

lim  <sN°  =  s(180°)  ■  11.4  cm. 

«  A  >0  ^  * 


This  is  equivalent  to  computing  the  undulation  correction  from  Stokes'  equation 
in  which  the  Stokes'  function  S(iji)  was  replaced  by  the  constant  $(180°). 


in  use. 

Of  ill 


However,  for  all  practical  applications  small  spherical  caps  are 
e.  In  our  example  the  truncation  effects  for  some  typical  values 


truncation  angle 


truncation  effect 


unmodified  kernel 

0.4  cm 
0.8  cm 
1. 1  cm 
1.5  cm 


modified  kernel 

0.2  cm 
0.4  cm 
0.6  cm 
0.9  cm 


Truncation  effect  on  the  zero-degree  undulation  term 
(formulae  35  and  36)  computed  for  the  zero-degree 
component  of  gravity  anomaly  Ag0  (eq.  34)  due  to  vertical 
datum  inconsistency  model  Ah3. 

Units  are  cm. 
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Figure  33.  Truncation  effect  (absolute  value)  on  the  zero-degree 
undulation  term  computed  according  to  formula  38,  using 
the  zero-degree  component  of  gravity  anomaly  Ag0  (eq.  34) 
due  to  vertical  datum  inconsistency  model  Ah3. 

Units  are  cm. 


We  see  that  the  modified  kernel  produces  much  smaller  truncation  effects 
than  the  unmodified  one.  In  fact  up  to  i^c=15°  the  effect  of  unmodified 
procedure  can  be  twice  the  effect  due  to  the  modified  one. 

We  can  also  evaluate  the  total  magnitude  of  the  effect  on  the 
zero-degree  undulation  according  to  formula  (31).  We  will  not  consider 


km 

2GR 


term,  since  6M  is  not  known.  The  second  term  in  equation  (31)  can 
be  evaluated  using  the  constants  of  GRS1980  used  in  previous  sections. 

G  =  9.797644  (m/s2) 

R  =  6371008.7714  m 

The  Ag0  effect  due  to  the  inconsistencies  of  vertical  datum  was  computed 
in  the  previous  chapter  as  0.57  x  10-2  mgal  (eq.  34).  Using  the  above 
constants  we  first  find: 

Ago  =  _1-85  cm 

Now  we  can  evaluate  the  joint  effect  of  AgQ  on  Nq  (eq.  31),  that 


4‘9o(1-q.o>  (38) 

for  different  values  of  the  truncation  radius  ^c.  Figure  33  shows 
this  effect  for  the  case  the  simple  Stokes'  kernel  is  used  exactly  as  in 
the  basic  eq.  (28),  and  for  the  case  a  modified  kernel  is  used  in  the 
integral  term  of  formula  (28). 

Now  we  will  evaluate  the  effect  of  datum  inconsistency,  which 
exists  in  a  given  set  of  terrestrial  gravity  anomalies  for  the  case 
the  geoid  undulations  are  computed  according  to  equation  (28)  using 
idealized  values  of  Laplace  harmonics  of  gravity  anomaly  field  implied 
by  potential  coefficients  (equation  28)  and  a  given  set  of  surface 
gravity  anomalies  in  a  cap  possibly  contaminated  by  vertical  datums 
inconsistency.  In  other  words  we  postulate  that  a  given  set  of  ter¬ 
restrial  gravity  anomalies  data  can  be  split  into  two  parts: 


Ag  =  Agj  +  6Ag  (39) 

where 

Ag j  -  idealized,  consistent  set  of  gravity  anomalies 

<5 Ag  -  the  effect  due  to  inconsistencies  in  vertical  datum. 

We  also  assume  that  the  idealized  set  of  terrestrial  data  ag,  is 
that  implied  by  the  potential  coefficients:  1 


oo 

A9t  =  l  Ag  (40) 

1  n=2  n 


where  agn  are  given  by  formula  (27). 

Substituting  for  Ag  in  equation  (28)  we  have 


N  =  No  +  4^G  ^  S(*}  {AgI  +  6Ag)  dCTc  +  C  Qln(  ,"c)  Agn  (41) 
o  n=2 

c 


or 


S(*>  49l  doc  +  S  "°  Qln(*c>  4dn  +  6N 
a  n=2 

c 

where 


6N  =  4?G  H  6Ag  dac 
a 

c 

Therefore,  we  can  rewrite  (42)  in  the  form: 


(42) 


(43) 


N  =  Nj  +  6N  (44) 

where 

Nj  -  idealized  set  of  geoid  undulation  implied  by  agj. 
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Now  we  will  evaluate  the  6N  term  (eq.  43)  which  can  be  attributed 
to  the  effect  of  datum  inconsistencies  implied  by  the  computational 
method  (28). 

Introducing  the  Molodenskii 1 s  truncation  coefficients  (eq.  24) 
we  can  rewrite  (43)  as: 


««  *  *TS  //  S<*>  “9  d°  -  35  l  0,„  649n  (45) 

a  n=2 


The  ordinary  Stokes'  integral  in  (45)  represents  the  total  effect  of 
datum  inconsistencies  in  terms  of  geoid  undulation  which  was  already 
evaluated  in  Chapter  6  (eq.  13)  and  is  given  in  Table  2  in  terms  of 
degree  variances  (we  will  use  data  in  the  third  column  of  Table  2  as 
pertaining  to  the  most  complete  third  model  of  datum  inconsistencies). 
The  quantities  6Ag  were  also  evaluated  in  Chapter  6  in  the  form  of 
degree  variances  (eq.  15).  They  simply  represent  the  spherical  har¬ 
monic  decomposition  of  the  effect  of  vertical  datum -inconsistencies 
described  in  terms  of  gravity  anomalies.  Now  we  can  rewrite  eq.  (45) 
using  the  explicit  spherical  harmonic  decomposition  of  the  integral 
term  in  formul^_(45)  and  substituting  for  6Agn  by  means  of  correction 
coefficients  6Cnm  and  6Snm  already  found  in  Chapter  6  (eq.  10): 

6N  =  1,  "  X  (‘E™»C°S  “  *  mX)  ?™»(cos  9)  * 

n-2  m-o 


-  I  Q-|n(0  I  (6CnmC0S  mA  +  6SnmS1n  Pnm(COS0)  (46) 

m  c  R  nm  nm  nm 

n-c  m-o 


or 

<5N  *  i  (R  l  (6CnmcosmA+  «Snmsin  mA)Pnm(cose  )(1-J1§-1-  Qln(x^c))  (47) 
n=2  m=0  u 


We  can  rewrite  (47)  as  follows: 

6N  =  j  Nn  (1  -  ¥  91n(*c)J  (48) 
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denotes  the  undulation  Laplace  spherical  harmonic  of  n-th  degree  due  to 
vertical  datum  inconsistencies.  Note  that  in  the  formula  (48)  the 
Laplace  harmonics  are  modified  by  the  cap-size  dependent  factor 


w> 


which  reflects  the  truncation  effect. 

Since  we  are  interested  mostly  in  the  magnitude  of  the  effect 
described  by  eq.  (48)  we  will  rewrite  this  equation  in  terms  of  degree 
variances  which  reflect  the  average  effect  due  to  the  specific  degree. 
Then,  the  cummulative  effect  up  to  degree  180  may  be  defined  as: 


Q  (*  ))2 

n  wc' ' 


(49) 


where 

xn  -  denotes  the  undulation  error  degree  variance  (Chapter  6,  eq.  13) 

The  numbers /x^j"  were  computed  in  Chapter  6  and  are  given  in  Table 
2.  We  will  use  the  last  column  in  Table  2  as  being  implied  by  our 
most  complete  model  of  datum  inconsistencies. 

Now,  equation  (49)  can  be  evaluated  for  different  values  of  cap- 
radius  <j ic  with  upper  summation  index  n  =180.  Figure  34  shows  this 
effect  plotted  for  the  two  cases:  Curve  1  was  constructed  exactly 
like  in  formula  (49)  using  unmodified  Molodenskii  ’ s  coefficients  QinU). 
Curve  2  was  constructed  using  modified  Molodenskii  coefficients  Q2nU) 
which  arise  when  the  modified  kernel  S(i|>)  -  S(ipc )  was  originally  used 
in  method  (28). 

From  Figure  34  we  see  that  maximum  effect  of  datum  inconsistencies 
occurs  for  the  truncation  angle  i^c=180°  (that  is  when  global  data  is  used), 
and  reaches  44.51  cm.  This  result  should  be  compared  with  44.7  cm 
cumulative  effect  on  geoid  undulation  computed  in  Chapter  6,  eq.  (14), 
which  is  based  on  the  summation  from  0  to  180  of  the  square  roots  of 
correction  coefficients  degree  variances  (eq.  10).  The  si ight  difference 
is  due  to  the  magnitude  of  the  zero  degree  term  which  is  not  considered 
by  formula  (49). 


ss 


Figure  34.  Cumulative  effect  on  geoid  undulation  (excluding  zero-degree 
term,  see  eq.  49)  due  to  vertical  datum  inconsistency  model 
Ah3  for  different  truncation  angles. 
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For  the  truncation  angle  y/c=0o  the  effect  due  to  eq.  (49)  is 
zero  since  no  terrestrial  gravity  data  s  actually  used.  For  spherical 
caps  having  radius  ipc  <  39°  the  effect  due  to  the  modified  Molodenskii 
kernel  is  considerably  smaller  than  the  effect  due  to  the  unmodified 
Stokes'  kernel. 

The  caps  used  in  most  practical  cases  have  radius  varying  from 
10  to  15  degrees  which  can  produce  the  effect  varying  from  6  to  11  cm 
for  the  case  of  modified  kernel  or  from  13  to  19  cm  for  the  case  of 
simple  Stokes'  kernel  (see  Figure  34). 

If  we  combine  the  cumulative  effect  on  undulation  due  to  harmonics 
having  degree  from  2  to  180  as  computed  from  formula  (49)  (Figure  34) 
with  the  zero  degree  effect  expressed  by  eq.  (38)  we  will  find  the 
total  effect  due  to  vertical  datum  inconsistencies  associated  with 
our  basic  method  for  computation  of  geoid  undulation  (see  eq.  28). 

Using  equations  (38)  and  (49)" the  magnitude  of  this  todal  effect  can 
be  expressed  as  follows: 


Figure  35  shows  this  effect  for  the  two  cases:  curve  1  is  related 
to  the  unmodified  Stokes'  kernel,  whereas  curve  2  is  related  to  the 
modified  Molodenskii ’s  kernel,  depending  which  one  is  used  in  our  basic 
method  of  geoid  undulation  computation  (formula  28).  From  the  figure 
we  note  that  for  the  truncation  angle  i^c=0°  the  total  effect  is  1.85  cm 
which  represents  the  constant: 


2G  Ago 


This  constant  is  independent  of  the  truncation  angle  ij»c.  It  was 
originally  introduced  by  the  formula  (22),  and  evaluated  for  the  zero 
degree  harmonic  component  of  the  error  in  gravity  anomaly  field  due 
to  the  inconsistency  in  the  vertical  datums  (see  formula  34).  For 
the  truncation  angle  ^c=180°,  that  is  for  the  global  data  set,  the 
total  effect  reaches  44.55  cm  for  the  ordinary  Stokes'  kernel  and 
45.52  cm  for  the  modified  one. 

For  the  usual  cap  size  ranging  from  10  to  15  degrees  the  total 
effect  shown  on  Figure  35  varies  from  13  to  19  cm  for  the  ordinary 
Stokes'  kernel  'and  from  7  to  11  cm  for  the  modified  one.  Since  in 
many  applications  there  is  a  need  for  the  determination  of  the  geoid 
to  the  accuracy  of  10  cm,  the  effects  computed  here  cannot  be  neglected. 
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Figure  35.  The  total  cumulative  effect  on  geoid  undulation  (including 
zero-degree  component,  see  eq.  50)  due  to  vertical  datum 
inconsistency  model  Ah3  for  different  truncation  angles. 


8.  SUMMARY  AND  CONCLUSIONS 


The  Sea  Surface  Topography  determines  the  deviation  of  the  Mean 
Sea  Level  from  the  geoid.  This  deviation  causes  different  vertical 
datums  existing  in  different  parts  of  the  world  to  refer  to  slightly 
different  zero-levels.  Those  zero-levels  define  different  equipotential 
surfaces  of  Earth's  gravitational  potential  or  surfaces  that  are  not 
even  equipotential  (due  to  overconstrained  adjustment  of  levelling 
net).  Therefore,  the  terrestrial  gravity  data  that  were  reduced  'to 
the  geoid'  are  in  fact  in  an  inconsistent  system.  These  inconsistencies 
in  the  global  vertical  datum  can  be  modeled  to  the  accuracy  the 
SST  is  known,  provided  sufficient  information  on  variety  of  existing 
vertical  networks  and  procedures  is  available. 

Today,  our  knowledge  on  the  SST  variation  is  far  from  complete. 

The  same  pertains  to  the  information  on  actual  vertical  datums  on  most 
of  the  land  areas  on  our  globe.  Therefore,  the  approach  presented 
in  this  work  can  be  classified  as  the  numerical  simulation  technique. 

We  constructed  three  different  approximations  of  the  hypothetic  vertical 
datum  inconsistencies  as  shown  on  Figures  17,  18  and  19.  In  order 
to  evaluate  the  propagated  error  in  potential  coefficients  due  to  the 
vertical  datum  inconsistencies  we  carried  out  the  spectral  analysis 
of  those  approximating  function  (formula  (10)).  The  correction  coef¬ 
ficients  produced  in  this  way  represent  the  errors  in  potential  coef¬ 
ficient  due  to  datum  inconsistency,  provided  geopotential  is  determined 
using  only  terrestrial  gravity  anomaly  data. 

Table  1  shows  that  global  datum  inconsistency  can  produce  a  sig¬ 
nificant  relative  error  in  potential  coefficients  that,  for  low  degree 
harmonics,  can  reach  almost  1  percent.  In  terms  of  geoid  undulation 
this  error  translates  to  the  average  of  almost  half  a  meter  at  first 
degree  harmonics.  The  cumulative  error  up  to  degree  180  (see  formula 
(14))  was  computed  as  44.67  cm,  which  means  it  can  easily  reach  the 
full  amplitude  of  the  SST  variation. 

Beyond  degree  60  the  e'fect  of  vertical  datum  inconsistency  is 
of  the  order  of  background  noise  and  can  be  neglected.  The  most  sig¬ 
nificant  contribution  is  in  the  low  degree  harmonics.  As  can  be  es¬ 
timated  from  Figure  30  the  low  order  correction  coefficients  within 
the  first  10  degrees  carry  as  much  as  92.6  percent  of  the  total  energy 
of  the  error-spectrum  (if  we  confine  our  study  up  to  degree  180). 

Beyond  degree  60  there  is  only  about  1  percent  of  the  total  energy 
in  the  error-spectrum. 

The  conclusion  is  that  by  supplying  the  low  degree  harmonics  of 
geopotential  (say  up  to  degree  10)  from  other  sources  than  the  ter¬ 
restrial  gravity  data  we  can  avoid  as  much  as  90  percent  of  the  error 
introduced  by  the  inconsistencies  in  vertical  datum.  From  that  point 
of  view  the  combination  of  satellite-derived  low  degree  coefficients 
with  the  high  degree  coefficients  derived  from  terrestrial  data  is 


highly  recommended.  Plots  similar  to  those  of  Figure  30  or  Figure  31 
can  provide  the  estimate  of  the  optimal  'splitting  degree',  so  that 
the  errors  over  the  whole  spectrum  are  minimized. 

The  internal  warp  in  the  North  Americal  vertical  datum  as  modeled 
in  Figure  15  seems  to  have  a  negligible  effect  on  the  determination 
of  spherical  harmonics  coefficients  of  gravitational  potential.  In 
terms  of  the  error  in  geoid  undulation,  this  effect  can  be  estimated 
to  be  of  the  order  of  millimeters. 

Finally,  the  effect  due  to  the  computational  procedures  such  as 
the  one  proposed  by  Rummel  and  Rapp  (1976)  (see  eq.  28)  is  of  a  great 
practical  importance.  It  was  shown  that  due  to  the  inconsistency  in 
the  world  vertical  datum  the  frequency  content  of  a  given  terrestrial 
gravity  anomaly  data  may  differ  from  the  frequency  content  of  an  ideal  ized 
anomalies  implied  by  a  consistent  vertical  reference.  This  discrepancy 
in  turn  will  propagate  on  geoid  undulation  as  determined  by  the  compu¬ 
tational  method  such  as  the  one  of  Rummel  and  Rapp  (eq.  28).  The  mag¬ 
nitude  of  that  effect  is  a  function  of  size  of  spherical  cap  used  in 
the  computations.  For  the  usual  cap  size  ranging  from  10  to  15  degrees 
in  spherical  radius  this  effect  ranges  from  7  cm  to  11  cm  or  from  13  cm 
to  19  cm  depending  on  the  modified  kernel  or  unmodified  one  was  used 
in  the  computations.  In  the  near  future  the  10  cm  accuracies  will 
be  required  for  the  determination  of  geoid.  We  conclude  that  due  to 
vertical  datums  inconsistencies  the  10  cm  level  is  actually  the  highest 
possible  resolution  that  could  be  achieved  using  the  terrestrial  gravity 
data.  One  cannot  go  beyond  that  limit  without  solving  the  vertical 
datum  problem. 
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